Spatial data analysis
Yazhuo Deng
January 06, 2019
- 0.0.1 Preliminaries about Spatial Data
- 0.0.2 Regression with Spatially Lagged Dependent Variables
- 0.0.3 The Spatial Error Model
- 0.0.4 Point-Referenced Data Models
- 0.0.5 Basics of Areal Data Models
- 0.0.6 Hierarchical Spatial Modeling (TODO)
- 0.0.7 Spatial Misalignment
- 0.0.8 Spatiotemporal Modeling
- 0.0.9 Spatial Point Pattern Analysis
0.0.1 Preliminaries about Spatial Data
#Spatial Data 在R上的图形展示
#https://zhuanlan.zhihu.com/p/27547697
#R code for Moran's I and Geary's C
#clean up
rm(list=ls())
#preferences
options(scipen=8)
#load required libraries
library(spdep) ## Loading required package: sp## Loading required package: Matrix## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(maps)
library(maptools)## Checking rgeos availability: TRUElibrary(foreign)
library(rgdal)## rgdal: version: 1.3-6, (SVN revision 773)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: FALSE
## Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
## Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.5/Resources/library/rgdal/proj
## Linking to sp version: 1.3-1#Load Data
imm.data.0<-read.dta("stateImmig0511.dta")
#fix(imm.data.0)
#names(imm.data.0)
#subset
imm.data<-subset(imm.data.0,State!="ALASKA" & State!="HAWAII")
#Estimate a Linear Model
mod.immigrant<-lm(immig0511~pubIdeolCCES, data=imm.data); summary(mod.immigrant)##
## Call:
## lm(formula = immig0511 ~ pubIdeolCCES, data = imm.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33009 -0.55823 0.04261 0.48461 1.38787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.93593 0.16676 5.613 0.0000010985 ***
## pubIdeolCCES 0.06738 0.01058 6.370 0.0000000806 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6671 on 46 degrees of freedom
## Multiple R-squared: 0.4687, Adjusted R-squared: 0.4571
## F-statistic: 40.58 on 1 and 46 DF, p-value: 0.00000008058#call data from "maps", create neighbor matrix
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0) #too cluttered, projection=???
IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs)
#eliminate Washington, DC
us.map<-us.map.0[-8]
#create neighbor list from state map
state.nb<-poly2nb(us.map)
state.coords<-coordinates(us.map)
plot(state.nb, coords=state.coords)
#state by state list, see how it compares to Melanie Wall's list: http://www.biostat.umn.edu/~brad/data/contig-lower48.dat
for(i in 1:48){print(state.nb[i])}## [[1]]
## [1] 8 9 22 40
##
## [[1]]
## [1] 4 5 26 29 42
##
## [[1]]
## [1] 16 22 23 34 40 41
##
## [[1]]
## [1] 2 26 35
##
## [[1]]
## [1] 2 14 25 29 34 42 48
##
## [[1]]
## [1] 19 30 37
##
## [[1]]
## [1] 18 28 36
##
## [[1]]
## [1] 1 9
##
## [[1]]
## [1] 1 8 31 38 40
##
## [[1]]
## [1] 24 26 35 42 45 48
##
## [[1]]
## [1] 12 13 15 23 47
##
## [[1]]
## [1] 11 15 20 33
##
## [[1]]
## [1] 11 21 23 25 39 47
##
## [[1]]
## [1] 5 23 25 34
##
## [[1]]
## [1] 11 12 23 33 40 44 46
##
## [[1]]
## [1] 3 22 41
##
## [[1]]
## [1] 27
##
## [[1]]
## [1] 7 36 44 46
##
## [[1]]
## [1] 6 27 30 37 43
##
## [[1]]
## [1] 12 33 47
##
## [[1]]
## [1] 13 32 39 47
##
## [[1]]
## [1] 1 3 16 40
##
## [[1]]
## [1] 3 11 13 14 15 25 34 40
##
## [[1]]
## [1] 10 32 39 48
##
## [[1]]
## [1] 5 13 14 23 39 48
##
## [[1]]
## [1] 2 4 10 35 42
##
## [[1]]
## [1] 17 19 43
##
## [[1]]
## [1] 7 30 36
##
## [[1]]
## [1] 2 5 34 41 42
##
## [[1]]
## [1] 6 19 28 36 43
##
## [[1]]
## [1] 9 38 40 44
##
## [[1]]
## [1] 21 24 39
##
## [[1]]
## [1] 12 15 20 36 46
##
## [[1]]
## [1] 3 5 14 23 29 41
##
## [[1]]
## [1] 4 10 26 45
##
## [[1]]
## [1] 7 18 28 30 33 46
##
## [[1]]
## [1] 6 19
##
## [[1]]
## [1] 9 31
##
## [[1]]
## [1] 13 21 24 25 32 48
##
## [[1]]
## [1] 1 3 9 15 22 23 31 44
##
## [[1]]
## [1] 3 16 29 34
##
## [[1]]
## [1] 2 5 10 26 29 48
##
## [[1]]
## [1] 19 27 30
##
## [[1]]
## [1] 15 18 31 40 46
##
## [[1]]
## [1] 10 35
##
## [[1]]
## [1] 15 18 33 36 44
##
## [[1]]
## [1] 11 13 20 21
##
## [[1]]
## [1] 5 10 24 25 39 42### MORAN'S I ###
#First, is the dependent variable spatially correlated on its own?
moran.test(imm.data$immig0511,nb2listw(state.nb,style='W'))##
## Moran I test under randomisation
##
## data: imm.data$immig0511
## weights: nb2listw(state.nb, style = "W")
##
## Moran I statistic standard deviate = 3.7071, p-value = 0.0001048
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.34159785 -0.02127660 0.00958191#Moran's I on residuals
moran.test(resid(mod.immigrant),nb2listw(state.nb,style='W'))##
## Moran I test under randomisation
##
## data: resid(mod.immigrant)
## weights: nb2listw(state.nb, style = "W")
##
## Moran I statistic standard deviate = 0.42017, p-value = 0.3372
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.019850914 -0.021276596 0.009581157# Moran's I test for residuals from a regression model
# based on weighted residuals
lm.morantest(mod.immigrant,nb2listw(state.nb,style='W'))##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = immig0511 ~ pubIdeolCCES, data = imm.data)
## weights: nb2listw(state.nb, style = "W")
##
## Moran I statistic standard deviate = 0.54009, p-value = 0.2946
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I Expectation Variance
## 0.019850914 -0.031418848 0.009011362### GEARY'S C ###
#Raw dependent variable
geary.test(imm.data$immig0511,nb2listw(state.nb,style='W'))##
## Geary C test under randomisation
##
## data: imm.data$immig0511
## weights: nb2listw(state.nb, style = "W")
##
## Geary C statistic standard deviate = 4.0725, p-value = 0.00002326
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.59254949 1.00000000 0.01001005#Geary's C on residuals
geary.test(resid(mod.immigrant),nb2listw(state.nb,style='W'))##
## Geary C test under randomisation
##
## data: resid(mod.immigrant)
## weights: nb2listw(state.nb, style = "W")
##
## Geary C statistic standard deviate = 0.77845, p-value = 0.2182
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.92211018 1.00000000 0.01001147#For homework:
#A tip for setting-up your neighbor matrix for a state using a map from the "maps" library
iowa<-map("county", "iowa", fill=TRUE, plot=TRUE, resolution=0) 
IDs <- sapply(strsplit(iowa$names, ","), function(x) x[2])
ia.map.0<-map2SpatialPolygons(iowa, IDs=IDs)
#Also for homework:
#Type ?cos for more information on trigonometric functions0.0.2 Regression with Spatially Lagged Dependent Variables
###WORKING EXAMPLE FROM CHAPTER 2 OF WARD & GLEDITSCH 2008###
###FUTURE USE OF THIS CODE REQUIRES CITATION OF THEIR ORIGINAL WORK###
#Load Data
# sldv is the data.frame
# mdd2 is the minimum distance data.frame
source("https://spia.uga.edu/faculty_pages/monogan/teaching/spatial/chapter1data.R")
#fix(sldv)
head(sldv)## polygon.number place.name tla country.name population kilometer2
## 1 3 Afghanistan AFG Afghanistan 17250390 641869.19
## 2 6 Albania ALB Albania 3416945 28754.50
## 3 4 Algeria ALG Algeria 27459230 2320972.00
## 4 9 Angola ANG Angola 11527260 1252421.00
## 5 11 Argentina ARG Argentina 33796870 2781013.00
## 6 7 Armenia ARM Armenia 3377228 29872.46
## idnumber gdp.2002 democracy
## 1 3 20000000000 -10
## 2 6 4840000000 5
## 3 4 55900000000 -3
## 4 9 11200000000 -3
## 5 11 102000000000 8
## 6 7 2370000000 5#fix(mdd2)
head(mdd2)## ida idb year mindist
## 114585 AFG CHN 2002 0
## 114586 AFG IND 2002 241
## 114587 AFG IRN 2002 0
## 114588 AFG KYR 2002 116
## 114589 AFG KZK 2002 331
## 114590 AFG OMA 2002 730# Create a neighbor minimum distance matrix
# from the dyadic minimum distance data
# (Note that the cshapes library provides a
# convienent way to create a matrix directly)
mddmat <- matrix(9999,ncol=dim(sldv)[1],nrow=dim(sldv)[1]) #We choose 9999km as the default because it's a big number that indicates "not a neighbor."
dimnames(mddmat) <- list(c(sldv$tla),c(sldv$tla))
for(i in 1:dim(mdd2)[1]){
mddmat[mdd2$ida[i],mdd2$idb[i]] <- mdd2$mindist[i] #super cool stuff: call text entries to set matrix cells
}
mddmat[1:12,1:12]## AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ALB 9999 9999 9999 9999 9999 9999 9999 554 9999 9999 9999 9999
## ALG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 763
## ANG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARM 9999 9999 9999 9999 9999 9999 9999 9999 0 9999 9999 9999
## AUL 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## AUS 9999 554 9999 9999 9999 9999 9999 9999 9999 9999 369 9999
## AZE 9999 9999 9999 9999 9999 0 9999 9999 9999 9999 9999 9999
## BAH 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## BEL 9999 9999 9999 9999 9999 9999 9999 369 9999 9999 9999 9999
## BEN 9999 9999 763 9999 9999 9999 9999 9999 9999 9999 9999 9999# Create a binary matrix with a 200 km threshold
m200mat <- mddmat
m200mat[m200mat<=200] <- 1
m200mat[m200mat>200] <- 0
m200mat[1:12,1:12]## AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG 0 0 0 0 0 0 0 0 0 0 0 0
## ALB 0 0 0 0 0 0 0 0 0 0 0 0
## ALG 0 0 0 0 0 0 0 0 0 0 0 0
## ANG 0 0 0 0 0 0 0 0 0 0 0 0
## ARG 0 0 0 0 0 0 0 0 0 0 0 0
## ARM 0 0 0 0 0 0 0 0 1 0 0 0
## AUL 0 0 0 0 0 0 0 0 0 0 0 0
## AUS 0 0 0 0 0 0 0 0 0 0 0 0
## AZE 0 0 0 0 0 1 0 0 0 0 0 0
## BAH 0 0 0 0 0 0 0 0 0 0 0 0
## BEL 0 0 0 0 0 0 0 0 0 0 0 0
## BEN 0 0 0 0 0 0 0 0 0 0 0 0# Create a row standardized weights matrix
wmat <- m200mat/apply(m200mat,1,sum)
wmat[1:12,1:12]## AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ALB 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ALG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ANG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ARG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ARM 0 0 0 0 0 0.0000000 0 0 0.2 0 0 0
## AUL 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## AUS 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## AZE 0 0 0 0 0 0.1428571 0 0 0.0 0 0 0
## BAH 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## BEL 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## BEN 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0# calculate the spatial lag of democracy
democracy.spatial.lag <- as.vector(wmat%*%sldv$democracy)
democracy.spatial.lag[1:12]## [1] -4.5714286 7.0000000 -0.4285714 2.2000000 8.6000000
## [6] 3.0000000 8.5000000 8.5454545 1.7142857 -10.0000000
## [11] 9.8000000 2.4000000# Create a weights list object (a list object is required for spdep)
cmat.lw <- spdep::mat2listw(m200mat, row.names=row.names(m200mat))
# Extract a neighbor list
cmat.nb <- cmat.lw$neighbours
### REALLY NOT RECOMMENDED: OLS WITH NO SPATIAL INFORMATION ###
stupid.model<-lm(democracy~log(gdp.2002/population), data=sldv); summary(stupid.model)##
## Call:
## lm(formula = democracy ~ log(gdp.2002/population), data = sldv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.945 -4.942 2.557 4.544 9.251
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -9.6874 2.4258 -3.994 0.000099996 ***
## log(gdp.2002/population) 1.6780 0.3128 5.364 0.000000289 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.285 on 156 degrees of freedom
## Multiple R-squared: 0.1557, Adjusted R-squared: 0.1503
## F-statistic: 28.77 on 1 and 156 DF, p-value: 0.0000002892moran.test(resid(stupid.model),nb2listw(cmat.nb))##
## Moran I test under randomisation
##
## data: resid(stupid.model)
## weights: nb2listw(cmat.nb)
##
## Moran I statistic standard deviate = 7.7664, p-value = 4.038e-15
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.404450410 -0.006369427 0.002798125### NOT RECOMMENDED: OLS WITH A SPATIAL LAG ###
sldv$dem.spat.lag<-democracy.spatial.lag
silly.model<-lm(democracy~log(gdp.2002/population)+dem.spat.lag,
data=sldv); summary(silly.model)##
## Call:
## lm(formula = democracy ~ log(gdp.2002/population) + dem.spat.lag,
## data = sldv)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.2447 -3.5407 0.6011 2.8480 11.7758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.97495 2.07116 -2.402 0.01749 *
## log(gdp.2002/population) 0.75927 0.27872 2.724 0.00719 **
## dem.spat.lag 0.76171 0.08802 8.653 6e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.177 on 155 degrees of freedom
## Multiple R-squared: 0.4307, Adjusted R-squared: 0.4234
## F-statistic: 58.64 on 2 and 155 DF, p-value: < 2.2e-16moran.test(resid(silly.model),nb2listw(cmat.nb))##
## Moran I test under randomisation
##
## data: resid(silly.model)
## weights: nb2listw(cmat.nb)
##
## Moran I statistic standard deviate = -3.1255, p-value = 0.9991
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.171459851 -0.006369427 0.002789921#But do not chase all that glitters... How are we doing Gauss-Markov wise? Poorly.
### BETTER CHOICE: MLE WITH A SPATIAL LAG ###
sldv.fit <- spdep::lagsarlm(democracy ~ log(gdp.2002/population),
data=sldv, nb2listw(cmat.nb),
method="eigen", quiet=FALSE); summary(sldv.fit)##
## Spatial lag model
## Jacobian calculated using neighbourhood matrix eigenvalues
## Computing eigenvalues ...
##
## rho: -0.4815483 function value: -551.2436
## rho: 0.0843528 function value: -508.153
## rho: 0.4340989 function value: -492.6041
## rho: 0.6502539 function value: -491.8775
## rho: 0.565321 function value: -491.0999
## rho: 0.559797 function value: -491.1005
## rho: 0.5631131 function value: -491.0994
## rho: 0.5631432 function value: -491.0994
## rho: 0.5631453 function value: -491.0994
## rho: 0.5631453 function value: -491.0994
## rho: 0.5631453 function value: -491.0994
## rho: 0.5631453 function value: -491.0994##
## Call:spdep::lagsarlm(formula = democracy ~ log(gdp.2002/population),
## data = sldv, listw = nb2listw(cmat.nb), method = "eigen",
## quiet = FALSE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.2640 -4.1620 1.2861 3.5655 11.1175
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.20342 2.08233 -2.9791 0.0028912
## log(gdp.2002/population) 0.99877 0.27826 3.5894 0.0003315
##
## Rho: 0.56315, LR test value: 45.035, p-value: 0.000000000019358
## Asymptotic standard error: 0.07576
## z-value: 7.4333, p-value: 1.0592e-13
## Wald statistic: 55.254, p-value: 1.0592e-13
##
## Log likelihood: -491.0994 for lag model
## ML residual variance (sigma squared): 27.161, (sigma: 5.2116)
## Number of observations: 158
## Number of parameters estimated: 4
## AIC: 990.2, (AIC for lm: 1033.2)
## LM test for residual autocorrelation
## test value: 2.1028, p-value: 0.14703spdep::moran.test(resid(sldv.fit),nb2listw(cmat.nb))##
## Moran I test under randomisation
##
## data: resid(sldv.fit)
## weights: nb2listw(cmat.nb)
##
## Moran I statistic standard deviate = -0.45541, p-value = 0.6756
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.030459980 -0.006369427 0.002798233### CALCULATING EQUILIBRIUM EFFECTS ###
# Code to calculate equilibrium effect of changes in GDP per capita
# Create vector to store the estimate for each states
ee.est <- rep(NA,dim(sldv)[1])
# Assign the country name labels
names(ee.est) <- sldv$tla
# Create a null vector to use in loop
svec <- rep(0,dim(sldv)[1])
# Create a N x N identity matrix
eye <- matrix(0,nrow=dim(sldv)[1],ncol=dim(sldv)[1])
diag(eye) <- 1
# Loop over 1:n states and store effect of change in
# each state i in ee.est[i]
for(i in 1:length(ee.est)){
cvec <- svec
cvec[i] <- 1
res <- solve(eye - 0.56315 * wmat) %*% cvec * 0.99877
ee.est[i] <- res[i]
}
## Results by country: How much would democracy increase in a country if log(gdp.2002/population) increased by 1 in that country? ##
ee.est## AFG ALB ALG ANG ARG ARM AUL AUS
## 1.067492 1.064252 1.104560 1.071895 1.141127 1.046587 1.121907 1.068377
## AZE BAH BEL BEN BFO BHM BHU BLR
## 1.054909 1.063698 1.089260 1.085159 1.093464 1.103532 1.077306 1.065284
## BLZ BNG BOL BOS BOT BRA BUI BUL
## 1.120205 1.084170 1.121271 1.057845 1.086653 1.160317 1.077169 1.060156
## CAM CAN CAO CDI CEN CHA CHL CHN
## 1.076181 1.071441 1.097271 1.089857 1.059433 1.061780 1.090393 1.104231
## COL CON COS CRO CUB CYP CZR DEN
## 1.122435 1.080874 1.114834 1.057845 1.151169 1.055224 1.043784 1.059586
## DJI DOM DRC DRV ECU EGY EQG ERI
## 1.094217 1.095957 1.086186 1.100169 1.074091 1.059802 1.081722 1.092907
## EST ETH FIN FJI FRN GAB GAM GFR
## 1.076870 1.097106 1.067623 1.130221 1.093156 1.089852 1.107729 1.087953
## GHA GNB GRC GRG GUA GUI GUY HAI
## 1.084570 1.107729 1.056691 1.049064 1.120205 1.154450 1.122284 1.225806
## HON HUN IND INS IRE IRN IRQ ISR
## 1.126266 1.063546 1.131937 1.151724 1.090164 1.080878 1.064917 1.073833
## ITA JAM JOR JPN KEN KUW KYR KZK
## 1.070141 1.101928 1.079516 1.108985 1.073287 1.043256 1.059848 1.070691
## LAO LAT LBR LEB LES LIB LIT LUX
## 1.084320 1.069949 1.120316 1.077060 1.060203 1.070449 1.064518 1.069460
## MAA MAC MAL MAW MEX MLD MLI MON
## 1.087258 1.060525 1.094080 1.066915 1.108328 1.034501 1.080262 1.028817
## MOR MYA MZM NAM NEP NIC NIG NIR
## 1.090795 1.071990 1.105845 1.086653 1.077306 1.139634 1.086168 1.078120
## NOR NTH OMA PAK PAN PAR PER PHI
## 1.074484 1.089260 1.064920 1.046347 1.129586 1.080619 1.119956 1.104765
## PNG POL POR PRK QAT ROK RUM RUS
## 1.243478 1.065041 1.090252 1.073413 1.092376 1.121669 1.068240 1.089609
## RWA SAF SAL SAU SEN SIE SLO SLV
## 1.067805 1.162269 1.114806 1.116874 1.118326 1.099129 1.055019 1.052143
## SOM SPN SRI SUD SUR SWA SWD SWZ
## 1.098215 1.153850 1.050053 1.074267 1.090706 1.067901 1.080689 1.039917
## SYR TAJ TAW TAZ THI TKM TOG TRI
## 1.081751 1.073590 1.087880 1.088571 1.094611 1.056263 1.085159 1.090907
## TUN TUR UAE UGA UKG UKR URU USA
## 1.050776 1.071617 1.076321 1.077890 1.152739 1.065221 1.064233 1.145732
## UZB VEN YEM YUG ZAM ZIM
## 1.066881 1.162103 1.101142 1.069267 1.092730 1.096897hist(ee.est)
## If log(gdp.2002/population) increased by 1 in Russia, what would happen to democracy in other states (observation 120)? ##
cvec <- rep(0,dim(sldv)[1])
cvec[120] <- 1
# Store estimates for impact of change in Russia in rus.est
eye <- matrix(0,nrow=dim(sldv)[1],ncol=dim(sldv)[1])
diag(eye) <- 1
rus.est <- solve(eye - 0.56315 * wmat) %*% cvec*0.99877
# Find ten highest values of rus.est vector
rus.est <- round(rus.est,3)
rus.est <- data.frame(sldv$tla,rus.est)
rus.est[rev(order(rus.est$rus.est)),][1:10,]## sldv.tla rus.est
## RUS RUS 1.090
## PRK PRK 0.242
## JPN JPN 0.241
## MON MON 0.240
## FIN FIN 0.222
## EST EST 0.212
## NOR NOR 0.202
## LIT LIT 0.198
## LAT LAT 0.197
## ARM ARM 0.177## What if China had a democracy score of 10? How would other states change in democracy? ##
# China is observation 32
cvec <- rep(0,dim(sldv)[1])
cvec[32] <- 10
# Store estimates of change in China in chn.est
chn.est <- c(cbind (0, 0, wmat%*%cvec) %*% c(summary(sldv.fit)$Coef[,1],summary(sldv.fit)$rho))
chn.est <- round(chn.est,3)
# Find all states where non-zero impact
chn.est <- data.frame(sldv$tla,chn.est)
chn.est.2<-subset(chn.est,chn.est>0)
chn.est.2[rev(order(chn.est.2$chn.est)),]## sldv.tla chn.est
## 139 TAW 1.877
## 116 PRK 1.877
## 96 MON 1.877
## 101 NEP 1.408
## 15 BHU 1.408
## 108 PAK 1.126
## 81 LAO 1.126
## 79 KYR 1.126
## 18 BNG 1.126
## 153 UZB 0.939
## 141 THI 0.939
## 98 MYA 0.939
## 138 TAJ 0.804
## 67 IND 0.804
## 44 DRV 0.804
## 1 AFG 0.804
## 80 KZK 0.704
## 120 RUS 0.282#### TIPS FOR HOMEWORK ####
# (1) #
#You have to load "chapter1data.R" again to get the "mdd2" information. We have not use for "sldv", though. After running "chapter1data.R", to avoid confusion, I'd type:
rm(sldv)
# (2) #
neighborhood.0<-read.csv("neighbor2002subset.csv")
#After you load the "neighbor2002subset.csv" dataset (which I named "neighborhood.0"), run the following line of code:
modeled<-!is.na(neighborhood.0$polity2) &
!is.na(neighborhood.0$lngdppercap) &
!is.na(neighborhood.0$physint)
#Create a second copy of your data subsetting on the conditions of "modeled".
# (3) #
#When defnining "mddmat" use code like this:
included<-c(as.character(neighborhood.0$tla))
mddmat <- matrix(9999,ncol=dim(neighborhood.0)[1],
nrow=dim(neighborhood.0)[1])
dimnames(mddmat) <- list(included,included)
for(i in 1:dim(mdd2)[1]){
if(is.element(mdd2$ida[i],included) & is.element(mdd2$idb[i],included)){
mddmat[mdd2$ida[i],mdd2$idb[i]] <- mdd2$mindist[i] #super cool stuff: call text entries to set matrix cells
}
}
#Here, "neighborhood" is my subsetted data frame. Convoluted code, I know, but I'll try to explain.0.0.3 The Spatial Error Model
#clean up
rm(list=ls())
#preferences
options(scipen=8)
#load required libraries
library(spdep)
library(maps)
library(foreign)
library(rgdal)
#Load Data
imm.data.0<-read.dta("stateImmig0511.dta")
#names(imm.data.0)
#subset
imm.data<-subset(imm.data.0,State!="ALASKA" & State!="HAWAII")
#call data from "maps", create neighbor matrix
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0)
IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs)
#eliminate Washington, DC
us.map<-us.map.0[-8]
#create neighbor list from state map
state.nb<-poly2nb(us.map)
state.coords<-coordinates(us.map)
plot(state.nb, coords=state.coords)
###estimate a simple linear model###
mod.immigrant<-lm(immig0511~pubIdeolCCES, data=imm.data); summary(mod.immigrant)##
## Call:
## lm(formula = immig0511 ~ pubIdeolCCES, data = imm.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33009 -0.55823 0.04261 0.48461 1.38787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.93593 0.16676 5.613 0.0000010985 ***
## pubIdeolCCES 0.06738 0.01058 6.370 0.0000000806 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6671 on 46 degrees of freedom
## Multiple R-squared: 0.4687, Adjusted R-squared: 0.4571
## F-statistic: 40.58 on 1 and 46 DF, p-value: 0.00000008058###Geary's C on residuals###
geary.test(resid(mod.immigrant),nb2listw(state.nb,style='C'))##
## Geary C test under randomisation
##
## data: resid(mod.immigrant)
## weights: nb2listw(state.nb, style = "C")
##
## Geary C statistic standard deviate = 1.5279, p-value = 0.06327
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic Expectation Variance
## 0.8346617 1.0000000 0.0117107###plot residuals###
X <- c(mod.immigrant$residuals[1:20],mod.immigrant$residuals[20],mod.immigrant$residuals[21:48])
n.col <- 10
br <- c(-999,quantile(X,c(.1,.2,.3,.4,.5,.6,.7,.8,.9)), 999)
#shading<-gray((n.col-1):0/(n.col-1))
shading<- colorRampPalette(c("red","white","blue"))(10)
#shading<- colorRampPalette(c("red","blue"))(10)
X.grp<-findInterval(X, vec=br, rightmost.closed = TRUE, all.inside = TRUE)
X.shad<-shading[X.grp]
states.to.map<-c("alabama","arizona","arkansas","california",
"colorado","connecticut","delaware","florida",
"georgia","idaho","illinois","indiana","iowa",
"kansas","kentucky","louisiana","maine",
"maryland","massachusetts:main","michigan:north",
"michigan:south","minnesota","mississippi",
"missouri","montana","nebraska","nevada",
"new hampshire","new jersey","new mexico",
"new york:main","north carolina:main","north dakota",
"ohio","oklahoma","oregon","pennsylvania",
"rhode island","south carolina","south dakota","tennessee",
"texas","utah","vermont","virginia:main",
"washington:main","west virginia","wisconsin","wyoming")
map(database="state", region=states.to.map, fill=T,
plot=T, interior=T, col=X.shad, exact=T)
#Fit spatial error model
spat.err.immig<-errorsarlm(immig0511~pubIdeolCCES, data=imm.data,
nb2listw(state.nb,style='C'), method="eigen",
quiet=FALSE); summary(spat.err.immig)##
## Spatial autoregressive error model
##
## Jacobian calculated using neighbourhood matrix eigenvalues
## Computing eigenvalues ...
##
## lambda: -0.6478654 function: -51.85541 Jacobian: -2.10781 SSE: 22.33468
## lambda: -0.08548168 function: -47.84376 Jacobian: -0.0385715 SSE: 20.59826
## lambda: 0.2620906 function: -47.83341 Jacobian: -0.4034806 SSE: 20.2787
## lambda: 0.09021095 function: -47.58365 Jacobian: -0.04491613 SSE: 20.37085
## lambda: 0.09006783 function: -47.58366 Jacobian: -0.04477178 SSE: 20.37097
## lambda: 0.09335054 function: -47.58356 Jacobian: -0.04814291 SSE: 20.36803
## lambda: 0.1578035 function: -47.61883 Jacobian: -0.1404927 SSE: 20.31965
## lambda: 0.0934995 function: -47.58356 Jacobian: -0.04829886 SSE: 20.3679
## lambda: 0.09349818 function: -47.58356 Jacobian: -0.04829749 SSE: 20.3679
## lambda: 0.09349813 function: -47.58356 Jacobian: -0.04829743 SSE: 20.3679
## lambda: 0.09349816 function: -47.58356 Jacobian: -0.04829747 SSE: 20.3679
## lambda: 0.09349868 function: -47.58356 Jacobian: -0.04829801 SSE: 20.3679
## lambda: 0.09349837 function: -47.58356 Jacobian: -0.04829769 SSE: 20.3679
## lambda: 0.09349826 function: -47.58356 Jacobian: -0.04829756 SSE: 20.3679
## lambda: 0.09349821 function: -47.58356 Jacobian: -0.04829752 SSE: 20.3679
## lambda: 0.09349819 function: -47.58356 Jacobian: -0.0482975 SSE: 20.3679
## lambda: 0.09349817 function: -47.58356 Jacobian: -0.04829748 SSE: 20.3679
## lambda: 0.09349818 function: -47.58356 Jacobian: -0.04829749 SSE: 20.3679##
## Call:errorsarlm(formula = immig0511 ~ pubIdeolCCES, data = imm.data,
## listw = nb2listw(state.nb, style = "C"), method = "eigen",
## quiet = FALSE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.33868 -0.51886 0.02905 0.47115 1.37927
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.922768 0.170728 5.4049 0.0000000648425
## pubIdeolCCES 0.065708 0.010713 6.1335 0.0000000008597
##
## Lambda: 0.093498, LR test value: 0.14407, p-value: 0.70427
## Asymptotic standard error: 0.2068
## z-value: 0.45213, p-value: 0.65118
## Wald statistic: 0.20442, p-value: 0.65118
##
## Log likelihood: -47.58356 for error model
## ML residual variance (sigma squared): 0.42433, (sigma: 0.65141)
## Number of observations: 48
## Number of parameters estimated: 4
## AIC: 103.17, (AIC for lm: 101.31)#Fit spatial lag model
spat.lag.immig<-lagsarlm(immig0511~pubIdeolCCES, data=imm.data,
nb2listw(state.nb,style='C'), method="eigen",
quiet=FALSE); summary(spat.lag.immig)##
## Spatial lag model
## Jacobian calculated using neighbourhood matrix eigenvalues
## Computing eigenvalues ...
##
## rho: -0.6478654 function value: -57.83021
## rho: -0.08548168 function value: -48.34439
## rho: 0.2620906 function value: -47.17023
## rho: 0.2022487 function value: -47.04967
## rho: 0.1784245 function value: -47.04129
## rho: 0.1814846 function value: -47.04113
## rho: 0.1813312 function value: -47.04112
## rho: 0.1813374 function value: -47.04112
## rho: 0.1813374 function value: -47.04112
## rho: 0.1813374 function value: -47.04112
## rho: 0.1813374 function value: -47.04112
## rho: 0.1813374 function value: -47.04112##
## Call:lagsarlm(formula = immig0511 ~ pubIdeolCCES, data = imm.data,
## listw = nb2listw(state.nb, style = "C"), method = "eigen",
## quiet = FALSE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.34985 -0.41520 0.08061 0.49837 1.40920
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.862126 0.171690 5.0214 0.00000051292
## pubIdeolCCES 0.061133 0.011192 5.4621 0.00000004707
##
## Rho: 0.18134, LR test value: 1.2289, p-value: 0.26761
## Asymptotic standard error: 0.15927
## z-value: 1.1386, p-value: 0.25489
## Wald statistic: 1.2963, p-value: 0.25489
##
## Log likelihood: -47.04112 for lag model
## ML residual variance (sigma squared): 0.41246, (sigma: 0.64223)
## Number of observations: 48
## Number of parameters estimated: 4
## AIC: 102.08, (AIC for lm: 101.31)
## LM test for residual autocorrelation
## test value: 0.54053, p-value: 0.462210.0.4 Point-Referenced Data Models
###Replication of Scallop example in R from chapter 2 of Banerjee, Carlin, & Gelfand 2004###
rm(list=ls())
#input data
myscallops <- read.table("myscallops.txt", header=T)
#load libraries
library(akima) #for drawing image/surface plots
library(geoR) #for variogram analyses## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.7-5.2.1 (built on 2016-05-02) is now loaded
## --------------------------------------------------------------library(scatterplot3d) #for drop line plot
#create an image plot with contours
int.scp <- akima::interp(x=myscallops$long,
y=myscallops$lat,
z=myscallops$lgcatch)
image(int.scp, xlim=range(myscallops$long),
ylim=range(myscallops$lat))
contour(int.scp, add=T)
#create a surface plot
persp(int.scp, xlim=range(myscallops$long),
ylim=range(myscallops$lat))
#drop line plot
scatterplot3d(x=myscallops$long, y=myscallops$lat, z=myscallops$lgcatch,
pch=20, highlight.3d=TRUE, type='h')
#define a "geodata" object
obj <- cbind(myscallops$long, myscallops$lat, myscallops$lgcatch)
scallops.geo <- geoR::as.geodata(obj, coords.col=1:2, data.col=3)
#create an empirical semivariogram
scallops.var <- geoR::variog(scallops.geo, estimator.type="classical")## variog: computing omnidirectional variogramscallops.var## $u
## [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
## [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
##
## $v
## [1] 3.048744 4.361375 4.520836 4.924096 5.697811 5.175957 4.888780
## [8] 4.345750 4.455269 4.224356 4.460676 5.338326 7.391131
##
## $n
## [1] 669 1538 1566 1534 1410 1267 1011 724 465 329 218 109 37
##
## $sd
## [1] 5.736655 6.292244 5.600070 6.288254 7.001918 6.329819 6.570810
## [8] 5.975687 5.815976 4.899078 5.451433 5.928083 7.090025
##
## $bins.lim
## [1] 0.000000000001 0.223081002026 0.446162004053 0.669243006079
## [5] 0.892324008105 1.115405010132 1.338486012158 1.561567014185
## [9] 1.784648016211 2.007729018237 2.230810020264 2.453891022290
## [13] 2.676972024316 2.900053026343
##
## $ind.bin
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##
## $var.mark
## [1] 4.720595
##
## $beta.ols
## [1] 3.482623
##
## $output.type
## [1] "bin"
##
## $max.dist
## [1] 2.900053
##
## $estimator.type
## [1] "classical"
##
## $n.data
## [1] 148
##
## $lambda
## [1] 1
##
## $trend
## [1] "cte"
##
## $pairs.min
## [1] 2
##
## $nugget.tolerance
## [1] 1e-12
##
## $direction
## [1] "omnidirectional"
##
## $tolerance
## [1] "none"
##
## $uvec
## [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
## [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
##
## $call
## geoR::variog(geodata = scallops.geo, estimator.type = "classical")
##
## attr(,"class")
## [1] "variogram"#create a robust empirical semivariogram
scallops.var.robust <- geoR::variog(scallops.geo, estimator.type="modulus")## variog: computing omnidirectional variogram#plot the two forms of the semivariogram
par(mfrow=c(2,1)) #align separate plots on the same space, two rows by one column
plot(scallops.var)
plot(scallops.var.robust)
#Fit a parametric variogram model (assume an exponential model) using weighted least squares
scallops.var.fit <- variofit(scallops.var.robust,
ini.cov.pars = c(6.0,1.0),
cov.model="exponential",
fix.nugget=FALSE, nugget=1.0)## variofit: covariance model used is exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optimscallops.var.fit## variofit: model parameters estimated by WLS (weighted least squares):
## covariance model is: exponential
## parameter estimates:
## tausq sigmasq phi
## 0.0000 5.1110 0.2065
## Practical Range with cor=0.05 for asymptotic range: 0.6184923
##
## variofit: minimised weighted sum of squares = 4017.547#Fit the exponential model using MLE
scallops.lik.fit <- likfit(scallops.geo, ini.cov.pars=c(6.0,1.0),
cov.model = "exponential", trend = "cte",
fix.nugget = FALSE, nugget = 1.0,
nospatial = TRUE, lik.method = "ML") ## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.#note: "lik.method" is a correction of a typo from the text
scallops.lik.fit## likfit: estimated model parameters:
## beta tausq sigmasq phi
## "2.3748" "0.0947" "5.7675" "0.2338"
## Practical Range with cor=0.05 for asymptotic range: 0.700551
##
## likfit: maximised log-likelihood = -285.9#CLASSICAL KRIGING
#First we must create locations to predict#
#Source: http://sudipto.bol.ucla.edu/Software/Book.R
u.loci <- rbind(c(-73.0,39.5),c(-72.5,40.0),c(-72,39.0))
#visualize where we're predicting
plot(y=myscallops$lat, x=myscallops$long, pch='.')
points(u.loci, pch='+', col='red')
scallops.classical <- geoR::krige.conv(geodata = scallops.geo,
locations = u.loci, borders = NULL,
krige = krige.control(cov.pars=c(5.7675,0.2338))) ## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood#cov.pars=c(sigmasq,phi)
scallops.classical$predict## [1] 5.675546 4.972033 2.067120#Bayesian Kriging Model--Estimation and Prediction
scallops.bayes2 <- geoR::krige.bayes(scallops.geo,
locations = u.loci, #defaults to locations="no"
borders = NULL,
model = model.control(trend.d = "cte",
cov.model = "exponential"),
prior = prior.control(beta.prior = "flat",
sigmasq.prior = "reciprocal",
tausq.rel.prior = "uniform",
tausq.rel.discrete = seq(from=0.0, to=1.0, by=0.01)))## Warning in geoR::krige.bayes(scallops.geo, locations = u.loci, borders =
## NULL, : .Random.seed not initialised. Creating it with runif(1)## krige.bayes: model with constant mean
## krige.bayes: computing the discrete posterior of phi/tausq.rel
## krige.bayes: argument `phi.discrete` not provided, using default values
## krige.bayes: computing the posterior probabilities.
## Number of parameter sets: 5050
## 1, 101, 201, 301, 401, 501, 601, 701, 801, 901, 1001, 1101, 1201, 1301, 1401, 1501, 1601, 1701, 1801, 1901, 2001, 2101, 2201, 2301, 2401, 2501, 2601, 2701, 2801, 2901, 3001, 3101, 3201, 3301, 3401, 3501, 3601, 3701, 3801, 3901, 4001, 4101, 4201, 4301, 4401, 4501, 4601, 4701, 4801, 4901, 5001,
##
## krige.bayes: sampling from posterior distribution
## krige.bayes: sample from the (joint) posterior of phi and tausq.rel
## [,1] [,2] [,3] [,4] [,5] [,6]
## phi 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148
## tausq.rel 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## frequency 18.0000000 7.0000000 5.0000000 8.0000000 9.0000000 4.0000000
## [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## phi 0.928017 1.044019 1.160021 1.276023 1.392025 1.62403 1.740032
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
## frequency 4.000000 3.000000 3.000000 4.000000 4.000000 3.00000 2.000000
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## phi 1.856034 1.972036 2.088038 2.20404 2.320042 2.436045 2.552047
## tausq.rel 0.000000 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000
## frequency 2.000000 4.000000 2.000000 3.00000 5.000000 1.000000 4.000000
## [,21] [,22] [,23] [,24] [,25] [,26] [,27]
## phi 2.668049 2.784051 2.900053 3.016055 3.132057 3.364062 3.480064
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 4.000000 2.000000 1.000000 3.000000 1.000000 1.000000 1.000000
## [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## phi 3.596066 3.712068 3.82807 3.944072 4.060074 4.176076 4.292078
## tausq.rel 0.000000 0.000000 0.00000 0.000000 0.000000 0.000000 0.000000
## frequency 2.000000 2.000000 4.00000 1.000000 5.000000 1.000000 2.000000
## [,35] [,36] [,37] [,38] [,39] [,40] [,41]
## phi 4.408081 4.524083 4.640085 4.872089 4.988091 5.104093 5.220095
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 1.000000 2.000000 3.000000 2.000000 2.000000 2.000000 1.000000
## [,42] [,43] [,44] [,45] [,46] [,47] [,48]
## phi 5.336098 5.4521 5.568102 5.684104 5.800106 0.2320042 0.3480064
## tausq.rel 0.000000 0.0000 0.000000 0.000000 0.000000 0.0100000 0.0100000
## frequency 2.000000 1.0000 1.000000 4.000000 1.000000 12.0000000 19.0000000
## [,49] [,50] [,51] [,52] [,53] [,54]
## phi 0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0100000 0.0100000 0.0100000 0.0100000 0.010000 0.010000
## frequency 5.0000000 9.0000000 8.0000000 3.0000000 6.000000 7.000000
## [,55] [,56] [,57] [,58] [,59] [,60] [,61]
## phi 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032 1.856034
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.01000 0.010000 0.010000
## frequency 2.000000 2.000000 2.000000 2.000000 7.00000 5.000000 7.000000
## [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## phi 1.972036 2.088038 2.20404 2.320042 2.436045 2.552047 2.668049
## tausq.rel 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000 0.010000
## frequency 5.000000 6.000000 6.00000 4.000000 3.000000 4.000000 1.000000
## [,69] [,70] [,71] [,72] [,73] [,74] [,75]
## phi 2.784051 2.900053 3.016055 3.132057 3.248059 3.364062 3.480064
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000
## frequency 4.000000 6.000000 1.000000 1.000000 3.000000 2.000000 5.000000
## [,76] [,77] [,78] [,79] [,80] [,81] [,82]
## phi 3.596066 3.712068 3.82807 4.060074 4.176076 4.292078 4.408081
## tausq.rel 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000 0.010000
## frequency 4.000000 2.000000 2.00000 2.000000 4.000000 2.000000 5.000000
## [,83] [,84] [,85] [,86] [,87] [,88] [,89]
## phi 4.524083 4.640085 4.872089 4.988091 5.220095 5.336098 5.4521
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.0100
## frequency 2.000000 3.000000 1.000000 1.000000 1.000000 5.000000 1.0000
## [,90] [,91] [,92] [,93] [,94] [,95]
## phi 5.568102 5.684104 5.800106 0.2320042 0.3480064 0.4640085
## tausq.rel 0.010000 0.010000 0.010000 0.0200000 0.0200000 0.0200000
## frequency 3.000000 1.000000 1.000000 16.0000000 11.0000000 16.0000000
## [,96] [,97] [,98] [,99] [,100] [,101]
## phi 0.5800106 0.6960127 0.8120148 0.928017 1.044019 1.160021
## tausq.rel 0.0200000 0.0200000 0.0200000 0.020000 0.020000 0.020000
## frequency 10.0000000 11.0000000 6.0000000 12.000000 3.000000 3.000000
## [,102] [,103] [,104] [,105] [,106] [,107] [,108]
## phi 1.276023 1.392025 1.508028 1.62403 1.740032 1.856034 1.972036
## tausq.rel 0.020000 0.020000 0.020000 0.02000 0.020000 0.020000 0.020000
## frequency 6.000000 5.000000 2.000000 2.00000 8.000000 5.000000 5.000000
## [,109] [,110] [,111] [,112] [,113] [,114] [,115]
## phi 2.088038 2.20404 2.320042 2.436045 2.552047 2.668049 2.900053
## tausq.rel 0.020000 0.02000 0.020000 0.020000 0.020000 0.020000 0.020000
## frequency 2.000000 3.00000 2.000000 5.000000 2.000000 2.000000 2.000000
## [,116] [,117] [,118] [,119] [,120] [,121]
## phi 3.016055 3.132057 3.596066 4.640085 0.1160021 0.2320042
## tausq.rel 0.020000 0.020000 0.020000 0.020000 0.0300000 0.0300000
## frequency 2.000000 1.000000 1.000000 1.000000 2.0000000 18.0000000
## [,122] [,123] [,124] [,125] [,126] [,127]
## phi 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148 0.928017
## tausq.rel 0.0300000 0.0300000 0.0300000 0.0300000 0.0300000 0.030000
## frequency 17.0000000 14.0000000 5.0000000 7.0000000 4.0000000 4.000000
## [,128] [,129] [,130] [,131] [,132] [,133] [,134]
## phi 1.044019 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032
## tausq.rel 0.030000 0.030000 0.030000 0.030000 0.030000 0.03000 0.030000
## frequency 3.000000 6.000000 2.000000 4.000000 1.000000 4.00000 3.000000
## [,135] [,136] [,137] [,138] [,139] [,140] [,141]
## phi 1.856034 1.972036 2.20404 2.552047 2.784051 3.364062 4.872089
## tausq.rel 0.030000 0.030000 0.03000 0.030000 0.030000 0.030000 0.030000
## frequency 1.000000 1.000000 1.00000 1.000000 2.000000 1.000000 1.000000
## [,142] [,143] [,144] [,145] [,146] [,147]
## phi 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148
## tausq.rel 0.0400000 0.0400000 0.0400000 0.0400000 0.0400000 0.0400000
## frequency 13.0000000 11.0000000 12.0000000 3.0000000 5.0000000 7.0000000
## [,148] [,149] [,150] [,151] [,152] [,153] [,154]
## phi 0.928017 1.044019 1.160021 1.276023 1.392025 1.62403 1.972036
## tausq.rel 0.040000 0.040000 0.040000 0.040000 0.040000 0.04000 0.040000
## frequency 6.000000 4.000000 1.000000 4.000000 1.000000 2.00000 1.000000
## [,155] [,156] [,157] [,158] [,159] [,160]
## phi 2.436045 2.784051 3.82807 0.2320042 0.3480064 0.4640085
## tausq.rel 0.040000 0.040000 0.04000 0.0500000 0.0500000 0.0500000
## frequency 1.000000 1.000000 1.00000 9.0000000 14.0000000 7.0000000
## [,161] [,162] [,163] [,164] [,165] [,166]
## phi 0.5800106 0.6960127 0.8120148 0.928017 1.044019 1.160021
## tausq.rel 0.0500000 0.0500000 0.0500000 0.050000 0.050000 0.050000
## frequency 4.0000000 9.0000000 3.0000000 3.000000 3.000000 4.000000
## [,167] [,168] [,169] [,170] [,171] [,172]
## phi 1.276023 1.392025 1.62403 1.740032 0.2320042 0.3480064
## tausq.rel 0.050000 0.050000 0.05000 0.050000 0.0600000 0.0600000
## frequency 1.000000 2.000000 1.00000 1.000000 14.0000000 9.0000000
## [,173] [,174] [,175] [,176] [,177] [,178]
## phi 0.4640085 0.5800106 0.8120148 0.928017 1.044019 1.276023
## tausq.rel 0.0600000 0.0600000 0.0600000 0.060000 0.060000 0.060000
## frequency 10.0000000 8.0000000 5.0000000 3.000000 4.000000 2.000000
## [,179] [,180] [,181] [,182] [,183] [,184]
## phi 0.1160021 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.0700000 0.0700000 0.0700000 0.0700000 0.0700000 0.0700000
## frequency 1.0000000 5.0000000 18.0000000 3.0000000 5.0000000 3.0000000
## [,185] [,186] [,187] [,188] [,189] [,190]
## phi 1.972036 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.070000 0.0800000 0.0800000 0.0800000 0.0800000 0.0800000
## frequency 1.000000 7.0000000 9.0000000 7.0000000 4.0000000 2.0000000
## [,191] [,192] [,193] [,194] [,195] [,196]
## phi 0.8120148 0.928017 1.044019 1.276023 1.740032 0.2320042
## tausq.rel 0.0800000 0.080000 0.080000 0.080000 0.080000 0.0900000
## frequency 1.0000000 3.000000 1.000000 1.000000 1.000000 9.0000000
## [,197] [,198] [,199] [,200] [,201] [,202]
## phi 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148 1.276023
## tausq.rel 0.0900000 0.0900000 0.0900000 0.0900000 0.0900000 0.090000
## frequency 8.0000000 7.0000000 6.0000000 2.0000000 2.0000000 1.000000
## [,203] [,204] [,205] [,206] [,207] [,208]
## phi 0.2320042 0.3480064 0.4640085 0.5800106 0.8120148 0.928017
## tausq.rel 0.1000000 0.1000000 0.1000000 0.1000000 0.1000000 0.100000
## frequency 2.0000000 6.0000000 6.0000000 2.0000000 1.0000000 1.000000
## [,209] [,210] [,211] [,212] [,213] [,214]
## phi 0.1160021 0.2320042 0.3480064 0.4640085 0.5800106 0.8120148
## tausq.rel 0.1100000 0.1100000 0.1100000 0.1100000 0.1100000 0.1100000
## frequency 1.0000000 3.0000000 4.0000000 2.0000000 3.0000000 1.0000000
## [,215] [,216] [,217] [,218] [,219] [,220]
## phi 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148
## tausq.rel 0.1200000 0.1200000 0.1200000 0.1200000 0.1200000 0.1200000
## frequency 7.0000000 5.0000000 1.0000000 4.0000000 1.0000000 1.0000000
## [,221] [,222] [,223] [,224] [,225] [,226]
## phi 0.2320042 0.3480064 0.4640085 0.5800106 0.2320042 0.3480064
## tausq.rel 0.1300000 0.1300000 0.1300000 0.1300000 0.1400000 0.1400000
## frequency 2.0000000 6.0000000 3.0000000 2.0000000 3.0000000 4.0000000
## [,227] [,228] [,229] [,230] [,231] [,232]
## phi 0.4640085 0.5800106 0.8120148 0.2320042 0.3480064 0.4640085
## tausq.rel 0.1400000 0.1400000 0.1400000 0.1500000 0.1500000 0.1500000
## frequency 1.0000000 1.0000000 1.0000000 2.0000000 4.0000000 1.0000000
## [,233] [,234] [,235] [,236] [,237] [,238]
## phi 0.2320042 0.3480064 0.6960127 0.2320042 0.3480064 0.4640085
## tausq.rel 0.1600000 0.1600000 0.1600000 0.1700000 0.1700000 0.1700000
## frequency 2.0000000 4.0000000 1.0000000 2.0000000 1.0000000 1.0000000
## [,239] [,240] [,241] [,242] [,243] [,244]
## phi 0.2320042 0.4640085 0.2320042 0.3480064 0.2320042 0.2320042
## tausq.rel 0.1800000 0.1800000 0.1900000 0.1900000 0.2000000 0.2100000
## frequency 6.0000000 4.0000000 2.0000000 4.0000000 1.0000000 1.0000000
## [,245] [,246] [,247] [,248] [,249] [,250]
## phi 0.5800106 0.2320042 0.2320042 0.4640085 0.3480064 0.4640085
## tausq.rel 0.2200000 0.2300000 0.2500000 0.2600000 0.2700000 0.2700000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 3.0000000 1.0000000
## [,251] [,252] [,253] [,254] [,255] [,256]
## phi 0.5800106 0.4640085 0.2320042 0.4640085 0.2320042 0.2320042
## tausq.rel 0.2700000 0.2900000 0.3100000 0.3600000 0.4100000 0.5300000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##
## krige.bayes: starting prediction at the provided locations
## krige.bayes: phi/tausq.rel samples for the predictive are same as for the posterior
## krige.bayes: computing moments of the predictive distribution
## krige.bayes: sampling from the predictive
## Number of parameter sets: 256
## 1, 11, 21, 31, 41, 51, 61, 71, 81, 91, 101, 111, 121, 131, 141, 151, 161, 171, 181, 191, 201, 211, 221, 231, 241, 251,
## krige.bayes: preparing summaries of the predictive distribution# Projected number of ads at three sites from Bayesian model
#Note: sites are set by: u.loci <- rbind(c(-73.0,39.5),c(-72.5,40.0),c(-72,39.0))
scallops.bayes2$predictive$mean## [1] 5.6620856 5.0163805 -0.3789692#Form the quantiles to interpret the model results
out <- scallops.bayes2$posterior
out <- out$sample
beta.qnt <- quantile(out$beta, c(0.50,0.025,0.975))
phi.qnt <-quantile(out$phi, c(0.50,0.025,0.975))
sigmasq.qnt <- quantile(out$sigmasq, c(0.50,0.025,0.975))
tausq.rel.qnt <- quantile(out$tausq.rel, c(0.50, 0.025, 0.975))
beta.qnt## 50% 2.5% 97.5%
## 1.821189 -6.790915 7.366864phi.qnt## 50% 2.5% 97.5%
## 0.5800106 0.2320042 4.9880912sigmasq.qnt## 50% 2.5% 97.5%
## 10.994736 4.226533 94.983472tausq.rel.qnt## 50% 2.5% 97.5%
## 0.03 0.00 0.18######################################################################
#R code for mapping locations of the 1974 pilot data
#clean up
rm(list=ls())
#load libraries
library(spatstat)## Loading required package: spatstat.data## Loading required package: nlme## Loading required package: rpart##
## spatstat 1.57-1 (nickname: 'Cartoon Physics')
## For an introduction to spatstat, type 'beginner'library(maptools)
library(rgdal)
library(maps)
library(foreign)
library(mgcv)## This is mgcv 1.8-26. For overview type 'help("mgcv-package")'.library(rgeos)## rgeos version: 0.4-2, (SVN revision 581)
## GEOS runtime version: 3.6.1-CAPI-1.10.1
## Linking to sp version: 1.3-1
## Polygon checking: TRUElibrary(spdep)
library(ggmap)## Loading required package: ggplot2checkCRSArgs <- function(uprojargs) {
.Call("checkCRSArgs", uprojargs, PACKAGE="rgdal")
}
#options
options(scipen=8)
#load the 1974 pilot data
pilot<-read.csv("pilot74.csv")
#visualize the locations
us.states<-map("state",fill=FALSE, plot=TRUE, resolution=0)
IDs <- us.states$names
us.states.1<-map(database="state",region=IDs, fill=TRUE,
plot=TRUE, interior=TRUE, exact=TRUE, resolution=0)
a.1<-maptools::map2SpatialPolygons(us.states.1, IDs=IDs,
proj4string=CRS("+proj=longlat"))
us.map.1<-spTransform(a.1,CRS("+init=esri:102003"))
###USEFUL WEBSITES###
# http://spatialreference.org/ref/esri/102003/
# https://www.nhgis.org/
plot(us.map.1,xlim=us.map.1@bbox[1,],ylim=us.map.1@bbox[2,])
axis(1);axis(2);box()
points(y=pilot$northings,x=pilot$eastings,col='red',pch='.')
0.0.5 Basics of Areal Data Models
#Estimating CAR models with world democracy data
#clean up
rm(list=ls())
#preferences
options(scipen=8)
#load required libraries
library(spdep)
library(maps)
library(foreign)
library(rgdal)
#Load Data
# sldv is the data.frame
# mdd2 is the minimum distance data.frame
source("https://spia.uga.edu/faculty_pages/monogan/teaching/spatial/chapter1data.R")
#fix(sldv)
#names(sldv)
#fix(mdd2)
#names(mdd2)
# Create a neighbor minimum distance matrix
# from the dyadic minimum distance data
# (Note that the cshapes library provides a
# convienent way to create a matrix directly)
mddmat <- matrix(9999,ncol=dim(sldv)[1],nrow=dim(sldv)[1])
#We choose 9999km as the default because it's a big number that indicates "not a neighbor."
dimnames(mddmat) <- list(c(sldv$tla),c(sldv$tla))
for(i in 1:dim(mdd2)[1]){
mddmat[mdd2$ida[i],mdd2$idb[i]] <- mdd2$mindist[i]
#super cool stuff: call text entries to set matrix cells
}
mddmat[1:12,1:12]## AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ALB 9999 9999 9999 9999 9999 9999 9999 554 9999 9999 9999 9999
## ALG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 763
## ANG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARG 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## ARM 9999 9999 9999 9999 9999 9999 9999 9999 0 9999 9999 9999
## AUL 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## AUS 9999 554 9999 9999 9999 9999 9999 9999 9999 9999 369 9999
## AZE 9999 9999 9999 9999 9999 0 9999 9999 9999 9999 9999 9999
## BAH 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999 9999
## BEL 9999 9999 9999 9999 9999 9999 9999 369 9999 9999 9999 9999
## BEN 9999 9999 763 9999 9999 9999 9999 9999 9999 9999 9999 9999# Create a binary matrix with a 200 km threshold
m200mat <- mddmat
m200mat[m200mat<=200] <- 1
m200mat[m200mat>200] <- 0
m200mat[1:12,1:12]## AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG 0 0 0 0 0 0 0 0 0 0 0 0
## ALB 0 0 0 0 0 0 0 0 0 0 0 0
## ALG 0 0 0 0 0 0 0 0 0 0 0 0
## ANG 0 0 0 0 0 0 0 0 0 0 0 0
## ARG 0 0 0 0 0 0 0 0 0 0 0 0
## ARM 0 0 0 0 0 0 0 0 1 0 0 0
## AUL 0 0 0 0 0 0 0 0 0 0 0 0
## AUS 0 0 0 0 0 0 0 0 0 0 0 0
## AZE 0 0 0 0 0 1 0 0 0 0 0 0
## BAH 0 0 0 0 0 0 0 0 0 0 0 0
## BEL 0 0 0 0 0 0 0 0 0 0 0 0
## BEN 0 0 0 0 0 0 0 0 0 0 0 0# Create a row standardized weights matrix
wmat <- m200mat/apply(m200mat,1,sum)
wmat[1:12,1:12]## AFG ALB ALG ANG ARG ARM AUL AUS AZE BAH BEL BEN
## AFG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ALB 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ALG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ANG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ARG 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## ARM 0 0 0 0 0 0.0000000 0 0 0.2 0 0 0
## AUL 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## AUS 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## AZE 0 0 0 0 0 0.1428571 0 0 0.0 0 0 0
## BAH 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## BEL 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0
## BEN 0 0 0 0 0 0.0000000 0 0 0.0 0 0 0# calculate the spatial lag of democracy
democracy.spatial.lag <- as.vector(wmat%*%sldv$democracy)
democracy.spatial.lag[1:12]## [1] -4.5714286 7.0000000 -0.4285714 2.2000000 8.6000000
## [6] 3.0000000 8.5000000 8.5454545 1.7142857 -10.0000000
## [11] 9.8000000 2.4000000# Create a weights list object (a list object is required for spdep)
cmat.lw <- mat2listw(m200mat, row.names=row.names(m200mat))
# Extract a neighbor list
cmat.nb <- cmat.lw$neighbours
#Fit Conditional Autoregressive Error Model
sldv.car<-spautolm(democracy ~ log(gdp.2002/population), data=sldv, listw=nb2listw(cmat.nb), family='CAR'); summary(sldv.car)## Warning in spautolm(democracy ~ log(gdp.2002/population), data = sldv,
## listw = nb2listw(cmat.nb), : Non-symmetric spatial weights in CAR model##
## Call:
## spautolm(formula = democracy ~ log(gdp.2002/population), data = sldv,
## listw = nb2listw(cmat.nb), family = "CAR")
##
## Residuals:
## Min 1Q Median 3Q Max
## -16.84735 -3.37506 0.44561 2.94816 12.81756
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.17337 3.32570 -1.5556 0.119810
## log(gdp.2002/population) 1.19528 0.38461 3.1078 0.001885
##
## Lambda: 0.91466 LR test value: 50.195 p-value: 0.0000000000013922
## Numerical Hessian standard error of lambda: 0.061333
##
## Log likelihood: -488.5194
## ML residual variance (sigma squared): 24.496, (sigma: 4.9493)
## Number of observations: 158
## Number of parameters estimated: 4
## AIC: 985.04#Fit spatial error model from Ward & Gleditsch
sldv.err<-errorsarlm(democracy ~ log(gdp.2002/population), data=sldv, listw=nb2listw(cmat.nb), method='eigen'); summary(sldv.err)##
## Call:
## errorsarlm(formula = democracy ~ log(gdp.2002/population), data = sldv,
## listw = nb2listw(cmat.nb), method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -14.9490 -4.2318 1.4473 3.4322 11.5587
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.48651 3.07263 -2.4365 0.0148296
## log(gdp.2002/population) 1.38711 0.37911 3.6589 0.0002533
##
## Lambda: 0.58188, LR test value: 44.18, p-value: 0.000000000029953
## Asymptotic standard error: 0.076529
## z-value: 7.6034, p-value: 2.8866e-14
## Wald statistic: 57.811, p-value: 2.8866e-14
##
## Log likelihood: -491.5267 for error model
## ML residual variance (sigma squared): 27.14, (sigma: 5.2096)
## Number of observations: 158
## Number of parameters estimated: 4
## AIC: 991.05, (AIC for lm: 1033.2)#####################################################################
#Areal models and diagnostic tests with immigrant policy data
#Load Data
imm.data.0<-read.dta("stateImmig0511.dta")
#subset
imm.data<-subset(imm.data.0,State!="ALASKA" & State!="HAWAII")
#call data from "maps", create neighbor matrix
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0)
IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs)
#eliminate Washington, DC
us.map<-us.map.0[-8]
#create neighbor list from state map
state.nb<-poly2nb(us.map)
#Basic linear model
mod.immigrant<-lm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+
demUnif+pcgsp1000+termLimits+changeForeign, data=imm.data); summary(mod.immigrant)##
## Call:
## lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) +
## repUnif + demUnif + pcgsp1000 + termLimits + changeForeign,
## data = imm.data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.3148 -0.3415 0.1175 0.4027 0.7816
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.217430 0.725108 -1.679 0.10096
## pubIdeolCCES 0.037442 0.012468 3.003 0.00459 **
## sqrt(squireProfess) 1.643635 0.766027 2.146 0.03802 *
## repUnif 0.038726 0.037425 1.035 0.30700
## demUnif 0.021630 0.039944 0.542 0.59116
## pcgsp1000 0.032330 0.011810 2.737 0.00919 **
## termLimits -0.244108 0.189977 -1.285 0.20621
## changeForeign -0.004057 0.001470 -2.760 0.00868 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5746 on 40 degrees of freedom
## Multiple R-squared: 0.6572, Adjusted R-squared: 0.5972
## F-statistic: 10.95 on 7 and 40 DF, p-value: 0.0000001241#tests on this linear model
tests<-lm.LMtests(mod.immigrant,
listw=nb2listw(state.nb, style='W'),test="all");tests##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
##
## LMerr = 1.6814, df = 1, p-value = 0.1947
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
##
## LMlag = 0.11809, df = 1, p-value = 0.7311
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
##
## RLMerr = 2.1739, df = 1, p-value = 0.1404
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
##
## RLMlag = 0.61066, df = 1, p-value = 0.4345
##
##
## Lagrange multiplier diagnostics for spatial dependence
##
## data:
## model: lm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess)
## + repUnif + demUnif + pcgsp1000 + termLimits + changeForeign, data
## = imm.data)
## weights: nb2listw(state.nb, style = "W")
##
## SARMA = 2.292, df = 2, p-value = 0.3179#Fit Conditional Autoregressive Error Model
hisp.car<-spautolm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
pcgsp1000+termLimits+changeForeign, data=imm.data,
listw=nb2listw(state.nb, style='C'), family='CAR')
summary(hisp.car)##
## Call: spautolm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) +
## repUnif + demUnif + pcgsp1000 + termLimits + changeForeign,
## data = imm.data, listw = nb2listw(state.nb, style = "C"),
## family = "CAR")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.25275 -0.36994 0.10529 0.34431 0.79470
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3881630 0.6414357 -2.1642 0.0304528
## pubIdeolCCES 0.0311108 0.0111077 2.8008 0.0050971
## sqrt(squireProfess) 2.0038527 0.6392588 3.1347 0.0017206
## repUnif 0.0396493 0.0339832 1.1667 0.2433185
## demUnif 0.0313643 0.0362122 0.8661 0.3864211
## pcgsp1000 0.0315273 0.0102941 3.0627 0.0021939
## termLimits -0.2702251 0.1643909 -1.6438 0.1002184
## changeForeign -0.0045352 0.0012220 -3.7113 0.0002062
##
## Lambda: -0.65412 LR test value: 1.4668 p-value: 0.22586
## Numerical Hessian standard error of lambda: 0.50925
##
## Log likelihood: -36.40544
## ML residual variance (sigma squared): 0.25519, (sigma: 0.50516)
## Number of observations: 48
## Number of parameters estimated: 10
## AIC: 92.811#Fit Simultaneous Autoregressive Error Model
hisp.sar<-spautolm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
pcgsp1000+termLimits+changeForeign, data=imm.data,
listw=nb2listw(state.nb, style='C'), family='SAR')
summary(hisp.sar)##
## Call: spautolm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) +
## repUnif + demUnif + pcgsp1000 + termLimits + changeForeign,
## data = imm.data, listw = nb2listw(state.nb, style = "C"),
## family = "SAR")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.283042 -0.276762 0.097122 0.362259 0.827288
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4337314 0.6344365 -2.2599 0.0238306
## pubIdeolCCES 0.0292941 0.0109771 2.6687 0.0076155
## sqrt(squireProfess) 2.1195678 0.6203145 3.4169 0.0006333
## repUnif 0.0393386 0.0337906 1.1642 0.2443480
## demUnif 0.0342368 0.0361459 0.9472 0.3435454
## pcgsp1000 0.0310867 0.0101159 3.0730 0.0021189
## termLimits -0.2778785 0.1608923 -1.7271 0.0841482
## changeForeign -0.0046582 0.0011823 -3.9397 0.00008157
##
## Lambda: -0.42474 LR test value: 1.8893 p-value: 0.16928
## Numerical Hessian standard error of lambda: 0.30305
##
## Log likelihood: -36.19418
## ML residual variance (sigma squared): 0.25469, (sigma: 0.50467)
## Number of observations: 48
## Number of parameters estimated: 10
## AIC: 92.388#Fit spatial error model from Ward & Gleditsch
hisp.err<-errorsarlm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
pcgsp1000+termLimits+changeForeign, data=imm.data,
listw=nb2listw(state.nb, style='C'), method='eigen')
summary(hisp.err)##
## Call:
## errorsarlm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) +
## repUnif + demUnif + pcgsp1000 + termLimits + changeForeign,
## data = imm.data, listw = nb2listw(state.nb, style = "C"),
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.283042 -0.276762 0.097122 0.362259 0.827288
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.4337314 0.6344365 -2.2599 0.0238306
## pubIdeolCCES 0.0292941 0.0109771 2.6687 0.0076155
## sqrt(squireProfess) 2.1195679 0.6203145 3.4169 0.0006333
## repUnif 0.0393386 0.0337906 1.1642 0.2443480
## demUnif 0.0342368 0.0361459 0.9472 0.3435454
## pcgsp1000 0.0310867 0.0101159 3.0730 0.0021189
## termLimits -0.2778785 0.1608923 -1.7271 0.0841482
## changeForeign -0.0046582 0.0011823 -3.9397 0.00008157
##
## Lambda: -0.42474, LR test value: 1.8893, p-value: 0.16928
## Asymptotic standard error: 0.22908
## z-value: -1.8541, p-value: 0.063718
## Wald statistic: 3.4379, p-value: 0.063718
##
## Log likelihood: -36.19418 for error model
## ML residual variance (sigma squared): 0.25469, (sigma: 0.50467)
## Number of observations: 48
## Number of parameters estimated: 10
## AIC: 92.388, (AIC for lm: 92.278)#Fit spatial lag model from Ward & Gleditsch
hisp.lag<-lagsarlm(immig0511~pubIdeolCCES+sqrt(squireProfess)+repUnif+demUnif+
pcgsp1000+termLimits+changeForeign, data=imm.data,
listw=nb2listw(state.nb, style='C'), method='eigen')
summary(hisp.lag)##
## Call:lagsarlm(formula = immig0511 ~ pubIdeolCCES + sqrt(squireProfess) +
## repUnif + demUnif + pcgsp1000 + termLimits + changeForeign,
## data = imm.data, listw = nb2listw(state.nb, style = "C"),
## method = "eigen")
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.30540 -0.37566 0.06357 0.36865 0.80318
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.2334590 0.6558271 -1.8808 0.0600034
## pubIdeolCCES 0.0394609 0.0118007 3.3440 0.0008259
## sqrt(squireProfess) 1.8376377 0.7026897 2.6151 0.0089189
## repUnif 0.0384534 0.0338370 1.1364 0.2557769
## demUnif 0.0232128 0.0361121 0.6428 0.5203560
## pcgsp1000 0.0326916 0.0106907 3.0580 0.0022285
## termLimits -0.2932435 0.1780584 -1.6469 0.0995795
## changeForeign -0.0045720 0.0013929 -3.2823 0.0010297
##
## Rho: -0.14321, LR test value: 0.7304, p-value: 0.39275
## Asymptotic standard error: 0.15683
## z-value: -0.91314, p-value: 0.36117
## Wald statistic: 0.83382, p-value: 0.36117
##
## Log likelihood: -36.77362 for lag model
## ML residual variance (sigma squared): 0.26979, (sigma: 0.51942)
## Number of observations: 48
## Number of parameters estimated: 10
## AIC: 93.547, (AIC for lm: 92.278)
## LM test for residual autocorrelation
## test value: 0.83283, p-value: 0.361460.0.6 Hierarchical Spatial Modeling (TODO)
https://aheblog.com/2016/12/07/geostatistical-modelling-with-r-and-stan/
Hierarchical Modeling and Analysis for Spatial Data, 2nd ed. (ISBN-13: 978-1-4398-1917-3), by S. Banerjee, B.P. Carlin and A.E. Gelfand, Boca Raton, FL: Chapman and Hall/CRC Press, 2015.
https://www.counterpointstat.com/hierarchical-modeling-and-analysis-for-spatial-data.html/
0.0.7 Spatial Misalignment
###Replication of Scallop example in R from chapter 2 of Banerjee, Carlin, & Gelfand 2004###
rm(list=ls())
#Create random observations within the block to krige
L <- 100
lat.min <- 39.5
lat.max <- 40.0
long.min <- -73.0
long.max <- -72.5
lat.vec <- runif(L,lat.min,lat.max)
long.vec <- runif(L,long.min,long.max)
#input data
myscallops <- read.table("myscallops.txt", header=T)
#myscallops <- read.table(file.choose(),header=T)
#load libraries
library(akima) #for drawing image/surface plots
library(geoR) #for variogram analyses
#create an image plot with contours #SKIP
int.scp <- interp(myscallops$long, myscallops$lat, myscallops$lgcatch)
image(int.scp, xlim=range(myscallops$long), ylim=range(myscallops$lat))
contour(int.scp, add=T)
polygon(x=c(-72.5,-72.5,-73,-73), y=c(40,39.5,39.5,40),border='cyan')
#What would you guess the average is?
#create a surface plot #SKIP
persp(int.scp, xlim=range(myscallops$long), ylim=range(myscallops$lat))
#define a "geodata" object
obj <- cbind(myscallops$long, myscallops$lat, myscallops$lgcatch)
scallops.geo <- as.geodata(obj, coords.col=1:2, data.col=3)
#create an empirical semivariogram
scallops.var <- variog(scallops.geo, estimator.type="classical")## variog: computing omnidirectional variogramscallops.var## $u
## [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
## [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
##
## $v
## [1] 3.048744 4.361375 4.520836 4.924096 5.697811 5.175957 4.888780
## [8] 4.345750 4.455269 4.224356 4.460676 5.338326 7.391131
##
## $n
## [1] 669 1538 1566 1534 1410 1267 1011 724 465 329 218 109 37
##
## $sd
## [1] 5.736655 6.292244 5.600070 6.288254 7.001918 6.329819 6.570810
## [8] 5.975687 5.815976 4.899078 5.451433 5.928083 7.090025
##
## $bins.lim
## [1] 0.000000000001 0.223081002026 0.446162004053 0.669243006079
## [5] 0.892324008105 1.115405010132 1.338486012158 1.561567014185
## [9] 1.784648016211 2.007729018237 2.230810020264 2.453891022290
## [13] 2.676972024316 2.900053026343
##
## $ind.bin
## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##
## $var.mark
## [1] 4.720595
##
## $beta.ols
## [1] 3.482623
##
## $output.type
## [1] "bin"
##
## $max.dist
## [1] 2.900053
##
## $estimator.type
## [1] "classical"
##
## $n.data
## [1] 148
##
## $lambda
## [1] 1
##
## $trend
## [1] "cte"
##
## $pairs.min
## [1] 2
##
## $nugget.tolerance
## [1] 1e-12
##
## $direction
## [1] "omnidirectional"
##
## $tolerance
## [1] "none"
##
## $uvec
## [1] 0.1115405 0.3346215 0.5577025 0.7807835 1.0038645 1.2269455 1.4500265
## [8] 1.6731075 1.8961885 2.1192695 2.3423505 2.5654315 2.7885125
##
## $call
## variog(geodata = scallops.geo, estimator.type = "classical")
##
## attr(,"class")
## [1] "variogram"#create a robust empirical semivariogram
scallops.var.robust <- variog(scallops.geo, estimator.type="modulus")## variog: computing omnidirectional variogram#plot the two forms of the semivariogram
par(mfrow=c(2,1)) #align separate plots on the same space, two rows by one column
plot(scallops.var)
plot(scallops.var.robust)
#Fit a parametric variogram model (assume an exponential model) using weighted least squares
scallops.var.fit <- variofit(scallops.var.robust, ini.cov.pars = c(1.0,2.0),
cov.model="exponential", fix.nugget=FALSE, nugget=1.0)## variofit: covariance model used is exponential
## variofit: weights used: npairs
## variofit: minimisation function used: optim## Warning in variofit(scallops.var.robust, ini.cov.pars = c(1, 2), cov.model
## = "exponential", : unreasonable initial value for sigmasq + nugget (too
## low)#Fit the exponential model using MLE
scallops.lik.fit <- likfit(scallops.geo, ini.cov.pars=c(1.0,2.0),
cov.model = "exponential", trend = "cte",
fix.nugget = FALSE, nugget = 1.0,
nospatial = TRUE, lik.method = "ML") ## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.#note: "lik.method" is a correction of a typo from the text
scallops.lik.fit## likfit: estimated model parameters:
## beta tausq sigmasq phi
## "2.3748" "0.0947" "5.7675" "0.2338"
## Practical Range with cor=0.05 for asymptotic range: 0.7005547
##
## likfit: maximised log-likelihood = -285.9#CLASSICAL KRIGING
#First we must create locations to predict#
#Source: http://www.biostat.umn.edu/~sudiptob/Software/Book.R
u.loci <- cbind(long.vec, lat.vec)
scallops.classical <- krige.conv(geodata = scallops.geo,
locations = u.loci,
borders = NULL,
krige = krige.control(cov.pars=c(5.7675,0.2338))) #cov.pars=c(sigmasq,phi)## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhoodscallops.classical$predict## [1] 5.6419341 1.7432806 6.9119948 2.7666652 6.5711829 4.1043280 0.7977967
## [8] 5.8868934 4.0250767 5.7365395 6.6456492 6.8274057 7.1243194 4.4532750
## [15] 6.1682690 6.9426476 3.9297373 3.8282980 6.1188330 6.3802473 6.7251522
## [22] 2.4312663 7.3685875 6.3593291 0.5356865 4.9771101 6.6788021 0.7907546
## [29] 5.6073281 5.9049236 4.1623659 4.4916631 6.3101120 5.4695400 6.6981826
## [36] 7.2099658 4.7542212 2.4198100 6.5279612 1.2411826 4.6001424 2.4825210
## [43] 5.6267846 3.5800146 3.1160066 4.0481120 3.1276008 7.6949181 6.5317052
## [50] 3.4223027 4.5079967 7.0502903 6.6481609 6.5888756 7.3429442 4.3821784
## [57] 2.3024721 4.8568017 5.2744735 5.5499560 1.8034411 5.2290168 7.1828775
## [64] 7.8271450 0.8225617 7.6170001 6.4242428 7.3011851 0.5264109 2.8854938
## [71] 2.9805537 4.3563635 6.8930906 7.2032837 6.3981411 7.0144640 2.4763012
## [78] 3.8205222 4.0216234 3.4520632 5.1477982 7.2604121 1.0833592 3.6994653
## [85] 3.6938089 5.3850278 4.4454305 4.7055334 5.8076492 3.6998101 7.0274229
## [92] 5.1036242 4.9551615 5.9113040 0.7688712 1.8666832 0.8744410 4.7890868
## [99] 3.7841463 7.4429925mean(scallops.classical$predict) #expected average scallop catch in the block## [1] 4.816904quantile(scallops.classical$predict, c(.5,.05,.95)) #median and 90% credible interval## 50% 5% 95%
## 5.0403671 0.8213234 7.3442264#Bayesian Kriging Model--Estimation and Prediction
scallops.bayes2 <- krige.bayes(scallops.geo,
locations = u.loci, #defaults to locations="no"
borders = NULL,
model = model.control(trend.d = "cte",
cov.model = "exponential"),
prior = prior.control(beta.prior = "flat",
sigmasq.prior = "reciprocal",
tausq.rel.prior = "uniform",
tausq.rel.discrete = seq(from=0.0, to=1.0, by=0.01)))## krige.bayes: model with constant mean
## krige.bayes: computing the discrete posterior of phi/tausq.rel
## krige.bayes: argument `phi.discrete` not provided, using default values
## krige.bayes: computing the posterior probabilities.
## Number of parameter sets: 5050
## 1, 101, 201, 301, 401, 501, 601, 701, 801, 901, 1001, 1101, 1201, 1301, 1401, 1501, 1601, 1701, 1801, 1901, 2001, 2101, 2201, 2301, 2401, 2501, 2601, 2701, 2801, 2901, 3001, 3101, 3201, 3301, 3401, 3501, 3601, 3701, 3801, 3901, 4001, 4101, 4201, 4301, 4401, 4501, 4601, 4701, 4801, 4901, 5001,
##
## krige.bayes: sampling from posterior distribution
## krige.bayes: sample from the (joint) posterior of phi and tausq.rel
## [,1] [,2] [,3] [,4] [,5] [,6]
## phi 0.1160021 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
## frequency 1.0000000 9.0000000 9.0000000 9.0000000 6.0000000 6.0000000
## [,7] [,8] [,9] [,10] [,11] [,12] [,13]
## phi 0.8120148 0.928017 1.044019 1.160021 1.276023 1.392025 1.508028
## tausq.rel 0.0000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 8.0000000 7.000000 2.000000 3.000000 4.000000 2.000000 1.000000
## [,14] [,15] [,16] [,17] [,18] [,19] [,20]
## phi 1.62403 1.740032 1.856034 1.972036 2.088038 2.20404 2.668049
## tausq.rel 0.00000 0.000000 0.000000 0.000000 0.000000 0.00000 0.000000
## frequency 4.00000 3.000000 1.000000 7.000000 2.000000 1.00000 2.000000
## [,21] [,22] [,23] [,24] [,25] [,26] [,27]
## phi 2.784051 3.016055 3.132057 3.248059 3.364062 3.480064 3.712068
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 4.000000 6.000000 2.000000 4.000000 2.000000 1.000000 3.000000
## [,28] [,29] [,30] [,31] [,32] [,33] [,34]
## phi 3.82807 3.944072 4.060074 4.176076 4.408081 4.524083 4.640085
## tausq.rel 0.00000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
## frequency 1.00000 1.000000 2.000000 3.000000 1.000000 5.000000 2.000000
## [,35] [,36] [,37] [,38] [,39] [,40] [,41]
## phi 4.756087 4.988091 5.104093 5.336098 5.4521 5.568102 0.2320042
## tausq.rel 0.000000 0.000000 0.000000 0.000000 0.0000 0.000000 0.0100000
## frequency 3.000000 1.000000 1.000000 2.000000 2.0000 2.000000 17.0000000
## [,42] [,43] [,44] [,45] [,46] [,47]
## phi 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148 0.928017
## tausq.rel 0.0100000 0.0100000 0.0100000 0.0100000 0.0100000 0.010000
## frequency 17.0000000 10.0000000 7.0000000 4.0000000 11.0000000 11.000000
## [,48] [,49] [,50] [,51] [,52] [,53] [,54]
## phi 1.044019 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.01000 0.010000
## frequency 5.000000 7.000000 6.000000 3.000000 4.000000 4.00000 4.000000
## [,55] [,56] [,57] [,58] [,59] [,60] [,61]
## phi 1.856034 1.972036 2.088038 2.20404 2.320042 2.436045 2.552047
## tausq.rel 0.010000 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000
## frequency 4.000000 4.000000 4.000000 7.00000 7.000000 5.000000 3.000000
## [,62] [,63] [,64] [,65] [,66] [,67] [,68]
## phi 2.668049 2.784051 2.900053 3.016055 3.132057 3.248059 3.364062
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000
## frequency 5.000000 3.000000 2.000000 5.000000 1.000000 2.000000 2.000000
## [,69] [,70] [,71] [,72] [,73] [,74] [,75]
## phi 3.480064 3.596066 3.82807 3.944072 4.060074 4.176076 4.292078
## tausq.rel 0.010000 0.010000 0.01000 0.010000 0.010000 0.010000 0.010000
## frequency 5.000000 2.000000 1.00000 3.000000 5.000000 5.000000 2.000000
## [,76] [,77] [,78] [,79] [,80] [,81] [,82]
## phi 4.408081 4.640085 4.756087 4.872089 4.988091 5.104093 5.336098
## tausq.rel 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000 0.010000
## frequency 2.000000 1.000000 2.000000 3.000000 4.000000 2.000000 1.000000
## [,83] [,84] [,85] [,86] [,87] [,88]
## phi 5.568102 5.684104 0.2320042 0.3480064 0.4640085 0.5800106
## tausq.rel 0.010000 0.010000 0.0200000 0.0200000 0.0200000 0.0200000
## frequency 1.000000 2.000000 17.0000000 10.0000000 14.0000000 15.0000000
## [,89] [,90] [,91] [,92] [,93] [,94]
## phi 0.6960127 0.8120148 0.928017 1.044019 1.160021 1.276023
## tausq.rel 0.0200000 0.0200000 0.020000 0.020000 0.020000 0.020000
## frequency 11.0000000 10.0000000 4.000000 3.000000 6.000000 7.000000
## [,95] [,96] [,97] [,98] [,99] [,100] [,101]
## phi 1.392025 1.508028 1.62403 1.740032 1.856034 1.972036 2.088038
## tausq.rel 0.020000 0.020000 0.02000 0.020000 0.020000 0.020000 0.020000
## frequency 2.000000 6.000000 4.00000 2.000000 5.000000 4.000000 5.000000
## [,102] [,103] [,104] [,105] [,106] [,107] [,108]
## phi 2.20404 2.320042 2.436045 2.552047 2.668049 2.784051 2.900053
## tausq.rel 0.02000 0.020000 0.020000 0.020000 0.020000 0.020000 0.020000
## frequency 2.00000 1.000000 4.000000 4.000000 2.000000 1.000000 1.000000
## [,109] [,110] [,111] [,112] [,113] [,114] [,115]
## phi 3.016055 3.248059 3.596066 3.712068 4.060074 4.408081 4.524083
## tausq.rel 0.020000 0.020000 0.020000 0.020000 0.020000 0.020000 0.020000
## frequency 1.000000 1.000000 1.000000 3.000000 1.000000 1.000000 1.000000
## [,116] [,117] [,118] [,119] [,120] [,121]
## phi 4.640085 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.020000 0.0300000 0.0300000 0.0300000 0.0300000 0.0300000
## frequency 1.000000 14.0000000 18.0000000 13.0000000 7.0000000 6.0000000
## [,122] [,123] [,124] [,125] [,126] [,127] [,128]
## phi 0.8120148 0.928017 1.044019 1.160021 1.276023 1.392025 1.508028
## tausq.rel 0.0300000 0.030000 0.030000 0.030000 0.030000 0.030000 0.030000
## frequency 10.0000000 6.000000 4.000000 5.000000 3.000000 2.000000 4.000000
## [,129] [,130] [,131] [,132] [,133] [,134] [,135]
## phi 1.62403 1.740032 1.856034 1.972036 2.20404 4.756087 0.2320042
## tausq.rel 0.03000 0.030000 0.030000 0.030000 0.03000 0.030000 0.0400000
## frequency 1.00000 2.000000 1.000000 4.000000 1.00000 1.000000 11.0000000
## [,136] [,137] [,138] [,139] [,140] [,141]
## phi 0.3480064 0.4640085 0.5800106 0.6960127 0.8120148 0.928017
## tausq.rel 0.0400000 0.0400000 0.0400000 0.0400000 0.0400000 0.040000
## frequency 11.0000000 6.0000000 9.0000000 7.0000000 6.0000000 9.000000
## [,142] [,143] [,144] [,145] [,146] [,147] [,148]
## phi 1.044019 1.160021 1.276023 1.392025 1.508028 1.62403 1.740032
## tausq.rel 0.040000 0.040000 0.040000 0.040000 0.040000 0.04000 0.040000
## frequency 5.000000 3.000000 5.000000 1.000000 4.000000 1.00000 1.000000
## [,149] [,150] [,151] [,152] [,153] [,154]
## phi 2.436045 3.016055 3.248059 3.712068 0.2320042 0.3480064
## tausq.rel 0.040000 0.040000 0.040000 0.040000 0.0500000 0.0500000
## frequency 1.000000 1.000000 1.000000 1.000000 10.0000000 8.0000000
## [,155] [,156] [,157] [,158] [,159] [,160]
## phi 0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0500000 0.0500000 0.0500000 0.0500000 0.050000 0.050000
## frequency 9.0000000 11.0000000 6.0000000 3.0000000 3.000000 2.000000
## [,161] [,162] [,163] [,164] [,165] [,166]
## phi 1.160021 1.276023 1.392025 1.740032 0.2320042 0.3480064
## tausq.rel 0.050000 0.050000 0.050000 0.050000 0.0600000 0.0600000
## frequency 1.000000 3.000000 2.000000 1.000000 8.0000000 13.0000000
## [,167] [,168] [,169] [,170] [,171] [,172]
## phi 0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0600000 0.0600000 0.0600000 0.0600000 0.060000 0.060000
## frequency 7.0000000 7.0000000 6.0000000 5.0000000 3.000000 1.000000
## [,173] [,174] [,175] [,176] [,177] [,178]
## phi 1.160021 1.276023 1.392025 1.856034 0.2320042 0.3480064
## tausq.rel 0.060000 0.060000 0.060000 0.060000 0.0700000 0.0700000
## frequency 1.000000 1.000000 2.000000 1.000000 11.0000000 12.0000000
## [,179] [,180] [,181] [,182] [,183] [,184]
## phi 0.4640085 0.5800106 0.6960127 0.8120148 0.928017 1.044019
## tausq.rel 0.0700000 0.0700000 0.0700000 0.0700000 0.070000 0.070000
## frequency 7.0000000 5.0000000 5.0000000 1.0000000 1.000000 3.000000
## [,185] [,186] [,187] [,188] [,189] [,190]
## phi 1.392025 1.508028 0.1160021 0.2320042 0.3480064 0.4640085
## tausq.rel 0.070000 0.070000 0.0800000 0.0800000 0.0800000 0.0800000
## frequency 1.000000 1.000000 1.0000000 5.0000000 8.0000000 8.0000000
## [,191] [,192] [,193] [,194] [,195] [,196]
## phi 0.5800106 0.6960127 0.8120148 1.276023 0.2320042 0.3480064
## tausq.rel 0.0800000 0.0800000 0.0800000 0.080000 0.0900000 0.0900000
## frequency 3.0000000 1.0000000 2.0000000 2.000000 1.0000000 9.0000000
## [,197] [,198] [,199] [,200] [,201] [,202]
## phi 0.4640085 0.5800106 0.8120148 0.928017 1.740032 0.2320042
## tausq.rel 0.0900000 0.0900000 0.0900000 0.090000 0.090000 0.1000000
## frequency 2.0000000 4.0000000 1.0000000 2.000000 1.000000 7.0000000
## [,203] [,204] [,205] [,206] [,207] [,208]
## phi 0.3480064 0.4640085 0.5800106 0.8120148 0.2320042 0.3480064
## tausq.rel 0.1000000 0.1000000 0.1000000 0.1000000 0.1100000 0.1100000
## frequency 4.0000000 1.0000000 2.0000000 1.0000000 4.0000000 5.0000000
## [,209] [,210] [,211] [,212] [,213] [,214]
## phi 0.4640085 0.5800106 0.6960127 0.8120148 0.928017 0.2320042
## tausq.rel 0.1100000 0.1100000 0.1100000 0.1100000 0.110000 0.1200000
## frequency 5.0000000 1.0000000 1.0000000 1.0000000 1.000000 4.0000000
## [,215] [,216] [,217] [,218] [,219] [,220]
## phi 0.3480064 0.4640085 0.5800106 0.8120148 0.928017 1.044019
## tausq.rel 0.1200000 0.1200000 0.1200000 0.1200000 0.120000 0.120000
## frequency 3.0000000 3.0000000 1.0000000 1.0000000 1.000000 1.000000
## [,221] [,222] [,223] [,224] [,225] [,226]
## phi 1.160021 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.120000 0.1300000 0.1300000 0.1300000 0.1300000 0.1300000
## frequency 1.000000 3.0000000 3.0000000 2.0000000 1.0000000 2.0000000
## [,227] [,228] [,229] [,230] [,231] [,232]
## phi 0.8120148 0.2320042 0.3480064 0.4640085 0.5800106 0.6960127
## tausq.rel 0.1300000 0.1400000 0.1400000 0.1400000 0.1400000 0.1400000
## frequency 1.0000000 2.0000000 3.0000000 4.0000000 1.0000000 1.0000000
## [,233] [,234] [,235] [,236] [,237] [,238]
## phi 0.8120148 0.928017 0.2320042 0.3480064 0.4640085 0.5800106
## tausq.rel 0.1400000 0.140000 0.1500000 0.1500000 0.1500000 0.1500000
## frequency 1.0000000 1.000000 5.0000000 1.0000000 3.0000000 1.0000000
## [,239] [,240] [,241] [,242] [,243] [,244]
## phi 0.6960127 0.2320042 0.3480064 0.4640085 0.2320042 0.3480064
## tausq.rel 0.1500000 0.1600000 0.1600000 0.1600000 0.1700000 0.1700000
## frequency 1.0000000 2.0000000 1.0000000 2.0000000 4.0000000 1.0000000
## [,245] [,246] [,247] [,248] [,249] [,250]
## phi 0.4640085 0.2320042 0.3480064 0.2320042 0.3480064 0.2320042
## tausq.rel 0.1700000 0.1800000 0.1800000 0.1900000 0.1900000 0.2000000
## frequency 1.0000000 2.0000000 2.0000000 2.0000000 1.0000000 1.0000000
## [,251] [,252] [,253] [,254] [,255] [,256]
## phi 0.3480064 0.2320042 0.3480064 0.2320042 0.3480064 0.3480064
## tausq.rel 0.2000000 0.2200000 0.2200000 0.2300000 0.2300000 0.2400000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
## [,257] [,258] [,259] [,260] [,261]
## phi 0.2320042 0.3480064 0.2320042 0.3480064 0.2320042
## tausq.rel 0.2800000 0.2800000 0.2900000 0.3200000 0.5400000
## frequency 1.0000000 1.0000000 1.0000000 1.0000000 1.0000000
##
## krige.bayes: starting prediction at the provided locations
## krige.bayes: phi/tausq.rel samples for the predictive are same as for the posterior
## krige.bayes: computing moments of the predictive distribution
## krige.bayes: sampling from the predictive
## Number of parameter sets: 261
## 1, 11, 21, 31, 41, 51, 61, 71, 81, 91, 101, 111, 121, 131, 141, 151, 161, 171, 181, 191, 201, 211, 221, 231, 241, 251, 261,
## krige.bayes: preparing summaries of the predictive distribution#TURN TO RESULTS PAGE
# Projected number of catches from Bayesian model
mean(scallops.bayes2$predictive$mean)## [1] 4.811959quantile(scallops.bayes2$predictive$mean, c(.5,.05,.95))## 50% 5% 95%
## 5.024128 1.248310 7.123376#Form the quantiles to interpret the model results
out <- scallops.bayes2$posterior
out <- out$sample
beta.qnt <- quantile(out$beta, c(0.50,0.025,0.975))
phi.qnt <-quantile(out$phi, c(0.50,0.025,0.975))
sigmasq.qnt <- quantile(out$sigmasq, c(0.50,0.025,0.975))
tausq.rel.qnt <- quantile(out$tausq.rel, c(0.50, 0.025, 0.975))
beta.qnt## 50% 2.5% 97.5%
## 1.766796 -7.892588 7.216852phi.qnt## 50% 2.5% 97.5%
## 0.6960127 0.2320042 4.7560870sigmasq.qnt## 50% 2.5% 97.5%
## 11.966236 4.258111 86.864834tausq.rel.qnt## 50% 2.5% 97.5%
## 0.03000 0.00000 0.160250.0.8 Spatiotemporal Modeling
###Replication of:
#Franzese and Hays. 2007. "Spatial Econometric Models of Cross-Sectional
#Interdependence in Political Science Panel and Time-Series-Cross-Section Data."
#Political Analysis 15:140-164.
###Usage of these data or code MUST cite the original paper.
#clean up
rm(list=ls())
#preferences
options(scipen=8)
#load required libraries
library(spdep)
library(maps)
library(foreign)
library(rgdal)
library(plm)## Loading required package: Formula#Load Data
fh<-read.csv('almMLire.csv',header=TRUE)
#Load Adjacency Matrix
adj<-as.matrix(read.csv('almWeights.csv',header=FALSE))
#Standardize Adjacency Matrix
noNeigh<-apply(adj,1,sum)
noNeigh[7]<-.01 #no neighbors for Greece, can't have 0 in the denominator
wmat <- adj/noNeigh
#Spatial Lag of the Outcome Variable, Using Kronecker Product
N<-15
T<-12
ident<-diag(1,nrow=T)
bigW<-(ident%x%wmat)
fh$SpatLag<-(ident%x%wmat)%*%fh$lnlmtue
#####REPLICATE TABLE 4#####
###First Model: Temporal Lag Only###
#dummied-out
mod.4.1.a<-lm(lnlmtue~lnlmtue_1+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+LLVOTE+
LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(cc)+as.factor(year),data=fh)
summary(mod.4.1.a)##
## Call:
## lm(formula = lnlmtue ~ lnlmtue_1 + DENSITY + DEIND + lngdp_pc +
## UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + OLDAGE +
## as.factor(cc) + as.factor(year), data = fh)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.19165 -0.14680 0.01142 0.12235 1.10187
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.3331830 3.5995002 -0.370 0.71165
## lnlmtue_1 0.5812769 0.0779578 7.456 0.00000000000804 ***
## DENSITY 0.0180919 0.0105857 1.709 0.08962 .
## DEIND 0.0619381 0.0338505 1.830 0.06938 .
## lngdp_pc -0.3283057 0.9457078 -0.347 0.72899
## UR -0.0892693 0.0334585 -2.668 0.00852 **
## TRADE 0.0007456 0.0085212 0.088 0.93039
## FDI -0.0228762 0.0249063 -0.918 0.35992
## LLVOTE -0.0130444 0.0153701 -0.849 0.39748
## LEFTC 0.0003791 0.0009957 0.381 0.70395
## TCDEMC -0.0025288 0.0035060 -0.721 0.47193
## GOVCON 0.0359730 0.0379986 0.947 0.34540
## OLDAGE -0.0106190 0.0585543 -0.181 0.85635
## as.factor(cc)2 -0.3317578 0.6096939 -0.544 0.58720
## as.factor(cc)3 -0.3808918 0.4829377 -0.789 0.43160
## as.factor(cc)4 -0.3877157 0.4384354 -0.884 0.37802
## as.factor(cc)5 0.5039088 0.4699241 1.072 0.28540
## as.factor(cc)6 0.7928225 0.2462236 3.220 0.00159 **
## as.factor(cc)7 0.1212794 0.4897753 0.248 0.80478
## as.factor(cc)8 0.0129930 0.5178942 0.025 0.98002
## as.factor(cc)9 0.0145534 0.5347751 0.027 0.97833
## as.factor(cc)10 -0.3887082 0.4065054 -0.956 0.34059
## as.factor(cc)11 -0.0886912 0.5295007 -0.167 0.86722
## as.factor(cc)12 0.3670486 0.5824373 0.630 0.52958
## as.factor(cc)13 -0.6791023 0.5635019 -1.205 0.23015
## as.factor(cc)14 0.4838704 0.4304817 1.124 0.26290
## as.factor(cc)15 -0.3977154 0.3402356 -1.169 0.24439
## as.factor(year)1988 0.1158841 0.1193023 0.971 0.33303
## as.factor(year)1989 0.3078653 0.1285753 2.394 0.01795 *
## as.factor(year)1990 0.1756465 0.1395881 1.258 0.21034
## as.factor(year)1991 -0.0388208 0.1502657 -0.258 0.79651
## as.factor(year)1992 0.0430456 0.1654094 0.260 0.79506
## as.factor(year)1993 0.0050429 0.1860829 0.027 0.97842
## as.factor(year)1994 0.0927372 0.2010415 0.461 0.64530
## as.factor(year)1995 0.1430783 0.2110666 0.678 0.49895
## as.factor(year)1996 0.1294835 0.2295124 0.564 0.57353
## as.factor(year)1997 0.1665405 0.2500523 0.666 0.50648
## as.factor(year)1998 0.2874073 0.2683866 1.071 0.28604
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3127 on 142 degrees of freedom
## Multiple R-squared: 0.8986, Adjusted R-squared: 0.8722
## F-statistic: 34.02 on 37 and 142 DF, p-value: < 2.2e-16#within estimator
mod.4.1.b<-plm(lnlmtue~lnlmtue_1+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+LLVOTE+
LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh,
index=c("cc","year"),model="within")
summary(mod.4.1.b)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = lnlmtue ~ lnlmtue_1 + DENSITY + DEIND + lngdp_pc +
## UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + OLDAGE +
## as.factor(year), data = fh, model = "within", index = c("cc",
## "year"))
##
## Balanced Panel: n = 15, T = 12, N = 180
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.191652 -0.146805 0.011423 0.122348 1.101871
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## lnlmtue_1 0.58127691 0.07795781 7.4563 0.000000000008038 ***
## DENSITY 0.01809190 0.01058574 1.7091 0.089620 .
## DEIND 0.06193812 0.03385047 1.8298 0.069384 .
## lngdp_pc -0.32830568 0.94570780 -0.3472 0.728990
## UR -0.08926930 0.03345849 -2.6681 0.008516 **
## TRADE 0.00074564 0.00852116 0.0875 0.930393
## FDI -0.02287622 0.02490634 -0.9185 0.359920
## LLVOTE -0.01304438 0.01537006 -0.8487 0.397484
## LEFTC 0.00037912 0.00099568 0.3808 0.703948
## TCDEMC -0.00252877 0.00350601 -0.7213 0.471932
## GOVCON 0.03597300 0.03799859 0.9467 0.345403
## OLDAGE -0.01061897 0.05855432 -0.1814 0.856350
## as.factor(year)1988 0.11588407 0.11930226 0.9713 0.333026
## as.factor(year)1989 0.30786531 0.12857529 2.3944 0.017951 *
## as.factor(year)1990 0.17564648 0.13958809 1.2583 0.210341
## as.factor(year)1991 -0.03882076 0.15026568 -0.2583 0.796513
## as.factor(year)1992 0.04304559 0.16540942 0.2602 0.795058
## as.factor(year)1993 0.00504291 0.18608285 0.0271 0.978418
## as.factor(year)1994 0.09273718 0.20104151 0.4613 0.645301
## as.factor(year)1995 0.14307828 0.21106656 0.6779 0.498949
## as.factor(year)1996 0.12948348 0.22951236 0.5642 0.573530
## as.factor(year)1997 0.16654050 0.25005231 0.6660 0.506477
## as.factor(year)1998 0.28740728 0.26838665 1.0709 0.286045
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.616
## Residual Sum of Squares: 13.881
## R-Squared: 0.58708
## Adj. R-Squared: 0.47949
## F-statistic: 8.77806 on 23 and 142 DF, p-value: < 2.22e-16###Second Model: OLS with a Spatial Lag###
mod.4.2<-plm(lnlmtue~lnlmtue_1+SpatLag+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+
LLVOTE+LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh,
index=c("cc","year"),model="within")
summary(mod.4.2)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = lnlmtue ~ lnlmtue_1 + SpatLag + DENSITY + DEIND +
## lngdp_pc + UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON +
## OLDAGE + as.factor(year), data = fh, model = "within", index = c("cc",
## "year"))
##
## Balanced Panel: n = 15, T = 12, N = 180
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.1573255 -0.1512225 0.0015748 0.1438301 1.0948248
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## lnlmtue_1 0.4354879852 0.0808440326 5.3868 0.0000002931 ***
## SpatLag -0.4376800180 0.1010789502 -4.3301 0.0000281051 ***
## DENSITY 0.0076330163 0.0102684684 0.7433 0.458510
## DEIND 0.0402701526 0.0323044963 1.2466 0.214618
## lngdp_pc -1.3073185100 0.9198428304 -1.4212 0.157455
## UR -0.0974263443 0.0316012054 -3.0830 0.002466 **
## TRADE 0.0000062976 0.0080356504 0.0008 0.999376
## FDI -0.0335940322 0.0236120494 -1.4227 0.157017
## LLVOTE 0.0043235810 0.0150359161 0.2876 0.774113
## LEFTC -0.0001629785 0.0009470531 -0.1721 0.863613
## TCDEMC -0.0028703711 0.0033064461 -0.8681 0.386807
## GOVCON 0.0224772985 0.0359607773 0.6251 0.532948
## OLDAGE 0.0516600127 0.0570484560 0.9055 0.366721
## as.factor(year)1988 0.2148278697 0.1147769654 1.8717 0.063320 .
## as.factor(year)1989 0.5480530853 0.1333104306 4.1111 0.0000665474 ***
## as.factor(year)1990 0.5016785281 0.1516218532 3.3087 0.001190 **
## as.factor(year)1991 0.2465674015 0.1562525590 1.5780 0.116805
## as.factor(year)1992 0.3011017340 0.1669491208 1.8036 0.073436 .
## as.factor(year)1993 0.2106993953 0.1817559899 1.1592 0.248316
## as.factor(year)1994 0.3030077992 0.1956656343 1.5486 0.123720
## as.factor(year)1995 0.3742281062 0.2060314733 1.8164 0.071439 .
## as.factor(year)1996 0.3771705533 0.2238195069 1.6852 0.094171 .
## as.factor(year)1997 0.4631084695 0.2454991847 1.8864 0.061298 .
## as.factor(year)1998 0.6604278507 0.2672999375 2.4707 0.014676 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.616
## Residual Sum of Squares: 12.252
## R-Squared: 0.63555
## Adj. R-Squared: 0.53733
## F-statistic: 10.2451 on 24 and 141 DF, p-value: < 2.22e-16###Third Model: Instrumenting the Spatial Lag in Two-Stage Least Squares###
#Spatially Lag the Xs
fh$l.DENSITY<-(ident%x%wmat)%*%fh$DENSITY
fh$l.DEIND<-(ident%x%wmat)%*%fh$DEIND
fh$l.lngdp_pc<-(ident%x%wmat)%*%fh$lngdp_pc
fh$l.UR<-(ident%x%wmat)%*%fh$UR
fh$l.TRADE<-(ident%x%wmat)%*%fh$TRADE
fh$l.FDI<-(ident%x%wmat)%*%fh$FDI
fh$l.LLVOTE<-(ident%x%wmat)%*%fh$LLVOTE
fh$l.LEFTC<-(ident%x%wmat)%*%fh$LEFTC
fh$l.TCDEMC<-(ident%x%wmat)%*%fh$TCDEMC
fh$l.GOVCON<-(ident%x%wmat)%*%fh$GOVCON
fh$l.OLDAGE<-(ident%x%wmat)%*%fh$OLDAGE
fh$l.lnlmtue_1<-(ident%x%wmat)%*%fh$lnlmtue_1
#Model the spatial lag of Y
#2SLS actually calls for modeling all the Xs! Hence, I'm getting somewhat different results.
stage.1<-lm(SpatLag~l.lnlmtue_1+l.DENSITY+l.DEIND+l.lngdp_pc+l.UR+l.TRADE+
l.FDI+l.LLVOTE+l.LEFTC+l.TCDEMC+l.GOVCON+l.OLDAGE, data=fh)#; summary(stage.1)#######
fh$iv.Lag<-stage.1$fitted.values
#Second Stage Equation
mod.4.3<-plm(lnlmtue~lnlmtue_1+iv.Lag+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+
LLVOTE+LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh,
index=c("cc","year"),model="within")
summary(mod.4.3)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = lnlmtue ~ lnlmtue_1 + iv.Lag + DENSITY + DEIND +
## lngdp_pc + UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON +
## OLDAGE + as.factor(year), data = fh, model = "within", index = c("cc",
## "year"))
##
## Balanced Panel: n = 15, T = 12, N = 180
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -1.2513435 -0.1325319 0.0022381 0.1336672 1.1433427
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## lnlmtue_1 0.49183370 0.08552898 5.7505 0.00000005301 ***
## iv.Lag -0.29385894 0.12417415 -2.3665 0.019317 *
## DENSITY 0.01133962 0.01080197 1.0498 0.295619
## DEIND 0.04422028 0.03414603 1.2950 0.197425
## lngdp_pc -0.74604335 0.94734235 -0.7875 0.432304
## UR -0.09255700 0.03295863 -2.8083 0.005688 **
## TRADE 0.00295506 0.00843820 0.3502 0.726711
## FDI -0.02923424 0.02465924 -1.1855 0.237803
## LLVOTE -0.00222941 0.01580224 -0.1411 0.888007
## LEFTC 0.00017133 0.00098386 0.1741 0.862006
## TCDEMC -0.00294343 0.00345501 -0.8519 0.395698
## GOVCON 0.03049498 0.03746921 0.8139 0.417093
## OLDAGE 0.02949127 0.06006905 0.4910 0.624221
## as.factor(year)1988 0.16224371 0.11903848 1.3630 0.175070
## as.factor(year)1989 0.40823609 0.13346051 3.0589 0.002660 **
## as.factor(year)1990 0.37013026 0.16008512 2.3121 0.022221 *
## as.factor(year)1991 0.19799317 0.17856377 1.1088 0.269400
## as.factor(year)1992 0.22534055 0.18009864 1.2512 0.212931
## as.factor(year)1993 0.17292861 0.19640025 0.8805 0.380092
## as.factor(year)1994 0.21630155 0.20463542 1.0570 0.292315
## as.factor(year)1995 0.28224374 0.21589190 1.3073 0.193226
## as.factor(year)1996 0.26475533 0.23300283 1.1363 0.257769
## as.factor(year)1997 0.31613859 0.25408695 1.2442 0.215485
## as.factor(year)1998 0.46706610 0.27483535 1.6994 0.091441 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.616
## Residual Sum of Squares: 13.35
## R-Squared: 0.60286
## Adj. R-Squared: 0.49583
## F-statistic: 8.91819 on 24 and 141 DF, p-value: < 2.22e-16###Fourth Model: Maximum Likelihood Estimator###
# For Comparison: Linear Model in MLE
linear.mle <- function(par, X, y, n) {
beta<-par[-length(par)]
sigma2<-par[length(par)]
loglikelihood <- -(n/2)*log(2*pi)-(n/2)*log(sigma2)-
(t(y-X%*%beta)%*%(y-X%*%beta)/(2*sigma2))
return(loglikelihood)
}
#subset data in matrix form
X.mod<-as.matrix(subset(fh,select=c(constant,lnlmtue_1,DENSITY,DEIND,lngdp_pc,
UR,TRADE,FDI,LLVOTE,LEFTC,TCDEMC,GOVCON,OLDAGE,
y88,y89,y90,y91,y92,y93,y94,y95,y96,y97,y98)))
# do a linear model with a time lag and time dummies
lin.mod <- optim(par=c(0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1), # starting value for beta and sigma
linear.mle, # the log-likelihood function
method="BFGS", # optimization method
hessian=TRUE, # return numerical Hessian
control=list(fnscale=-1), # maximize function instead of minimize
y=fh$lnlmtue, X=X.mod, n=180) # the data## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced# standard error of estimator
OI<-solve(-1*lin.mod$hessian)
se<-sqrt(diag(OI))
lin.mod$par; se## [1] 1.53859740656 0.75711795809 0.00075904866 0.00415123634
## [5] -0.01824582781 -0.02507694984 -0.00003919987 -0.01409629595
## [9] 0.00412294624 0.00095393272 0.00096184546 0.02718519240
## [13] -0.03505748214 0.09555748237 0.27821408818 0.11797066782
## [17] -0.09462528441 0.00215996393 -0.06591424593 0.05600699917
## [21] 0.13280965776 0.08864408076 0.14559411510 0.26519338209
## [25] 0.08645490262## [1] 0.7790569614 0.0489823918 0.0019540951 0.0135840395 0.1859839301
## [6] 0.0111310886 0.0016638735 0.0190482414 0.0098078452 0.0006761811
## [11] 0.0012633199 0.0112462704 0.0222351484 0.1079039004 0.1098943158
## [16] 0.1126536122 0.1126779789 0.1116047453 0.1138235388 0.1146197103
## [21] 0.1154613847 0.1186928176 0.1219921841 0.1395409277 0.0091100291#compare OLS to MLE
lin.mod.ols<-lm(lnlmtue~lnlmtue_1+DENSITY+DEIND+lngdp_pc+UR+TRADE+FDI+
LLVOTE+LEFTC+TCDEMC+GOVCON+OLDAGE+as.factor(year),data=fh)
summary(lin.mod.ols)##
## Call:
## lm(formula = lnlmtue ~ lnlmtue_1 + DENSITY + DEIND + lngdp_pc +
## UR + TRADE + FDI + LLVOTE + LEFTC + TCDEMC + GOVCON + OLDAGE +
## as.factor(year), data = fh)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.27217 -0.14718 0.01391 0.14856 1.19934
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.53815566 0.83681661 1.838 0.0679 .
## lnlmtue_1 0.75715953 0.05261397 14.391 <2e-16 ***
## DENSITY 0.00075865 0.00209897 0.361 0.7183
## DEIND 0.00415689 0.01459117 0.285 0.7761
## lngdp_pc -0.01822845 0.19977287 -0.091 0.9274
## UR -0.02508193 0.01195635 -2.098 0.0375 *
## TRADE -0.00003922 0.00178723 -0.022 0.9825
## FDI -0.01409620 0.02046049 -0.689 0.4919
## LLVOTE 0.00412575 0.01053500 0.392 0.6959
## LEFTC 0.00095433 0.00072631 1.314 0.1908
## TCDEMC 0.00096110 0.00135698 0.708 0.4798
## GOVCON 0.02717781 0.01208007 2.250 0.0259 *
## OLDAGE -0.03506607 0.02388367 -1.468 0.1441
## as.factor(year)1988 0.09556586 0.11590395 0.825 0.4109
## as.factor(year)1989 0.27819170 0.11804193 2.357 0.0197 *
## as.factor(year)1990 0.11791994 0.12100580 0.974 0.3313
## as.factor(year)1991 -0.09464141 0.12103198 -0.782 0.4354
## as.factor(year)1992 0.00218939 0.11987917 0.018 0.9855
## as.factor(year)1993 -0.06588254 0.12226247 -0.539 0.5908
## as.factor(year)1994 0.05602039 0.12311767 0.455 0.6497
## as.factor(year)1995 0.13276116 0.12402175 1.070 0.2861
## as.factor(year)1996 0.08858347 0.12749276 0.695 0.4882
## as.factor(year)1997 0.14550139 0.13103674 1.110 0.2685
## as.factor(year)1998 0.26512769 0.14988656 1.769 0.0789 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.3158 on 156 degrees of freedom
## Multiple R-squared: 0.8863, Adjusted R-squared: 0.8696
## F-statistic: 52.89 on 23 and 156 DF, p-value: < 2.2e-16## Log-Likelihood Fuction for Spatial and Time Lag##
spatiotemporal <- function(par, X, y, W, T, N) {
n<-N*T
rho<-par[1]
sigma2<-par[2]
beta<-par[3:length(par)]
ident<-diag(1,nrow=T)
nIdent<-diag(1,nrow=N)
bigIdent<-ident%x%nIdent
bigW<-(ident%x%W)
first<-prod(1-rho*eigen(bigW)$values) #Similar values for bigW or W. NaNs produced, though.
# first<-det(bigIdent-rho*bigW)
loglikelihood <- log(first)-(n/2)*log(2*pi)-(n/2)*log(sigma2)-
((t(y-rho*bigW%*%y-X%*%beta)%*%(y-rho*bigW%*%y-X%*%beta))/(2*sigma2))
return(loglikelihood)
}
#prep data for simple model
inds<-model.matrix(~factor(fh$cc)-1)
X.mod.b.0<-as.matrix(subset(fh,select=c(constant,lnlmtue_1,y88,y89,y90,
y91,y92,y93,y94,y95,y96,y97,y98)))
X.mod.b<-cbind(X.mod.b.0,inds)
#estimate simple model
spatiotemp.mod <- optim(par=c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0), # starting value for beta and sigma
spatiotemporal, # the log-likelihood function
method="BFGS", # optimization method
hessian=TRUE, # return numerical Hessian
control=list(fnscale=-1), # maximize function instead of minimize
y=fh$lnlmtue, X=X.mod.b,W=wmat, N=15, T=12) ## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs producedspatiotemp.mod$par## [1] -0.25897386 0.07752326 3.88603226 0.65694251 0.13879425
## [6] 0.39534813 0.28949429 0.03774227 0.08594637 -0.06105515
## [11] 0.03530374 0.11915186 0.07396601 0.13881230 0.26196905
## [16] 0.41322640 0.33317484 1.03345956 0.68676933 0.30464497
## [21] 0.57479921 -1.87122859 -0.03187805 0.55175184 0.74181189
## [26] -0.07687160 -0.11597257 0.90109413 0.35857262 0.08267825# standard error of estimator
OI<-solve(-1*spatiotemp.mod$hessian)
se<-sqrt(diag(OI))## Warning in sqrt(diag(OI)): NaNs producedspatiotemp.mod$par; se## [1] -0.25897386 0.07752326 3.88603226 0.65694251 0.13879425
## [6] 0.39534813 0.28949429 0.03774227 0.08594637 -0.06105515
## [11] 0.03530374 0.11915186 0.07396601 0.13881230 0.26196905
## [16] 0.41322640 0.33317484 1.03345956 0.68676933 0.30464497
## [21] 0.57479921 -1.87122859 -0.03187805 0.55175184 0.74181189
## [26] -0.07687160 -0.11597257 0.90109413 0.35857262 0.08267825## [1] 0.067119427 0.008221808 NaN 0.055020379 0.102583244
## [6] 0.107165680 0.113019875 0.110350365 0.106780185 0.104265401
## [11] 0.102911500 0.102835229 0.103713089 0.104187654 0.105700740
## [16] NaN NaN NaN NaN NaN
## [21] NaN NaN NaN NaN NaN
## [26] NaN NaN NaN NaN NaN##prep data for fuller model
inds<-model.matrix(~factor(fh$cc)-1)
X.mod<-as.matrix(subset(fh,select=c(constant,lnlmtue_1,DENSITY,DEIND,
lngdp_pc,UR,TRADE,FDI,LLVOTE,LEFTC,TCDEMC,GOVCON,
OLDAGE,y88,y89,y90,y91,y92,y93,y94,y95,y96,y97,y98)))
X.mod.a<-cbind(X.mod,inds)
#fit simple model
spatiotemp.full <- optim(par=c(0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0), # starting value for beta and sigma
spatiotemporal, # the log-likelihood function
method="BFGS", # optimization method
hessian=TRUE, # return numerical Hessian
control=list(fnscale=-1), # maximize function instead of minimize
y=fh$lnlmtue, X=X.mod.a,W=wmat, N=15, T=12) ## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs produced
## Warning in log(sigma2): NaNs producedspatiotemp.full$par## [1] -0.28958562211 0.06910988208 3.90211984452 0.48481701150
## [5] 0.01117189971 0.04760042241 -0.97604626466 -0.09466563444
## [9] 0.00025650397 -0.02996796219 -0.00155310619 0.00002044894
## [13] -0.00275479490 0.02704489833 0.03058603721 0.18134892304
## [17] 0.46678520401 0.39136404912 0.15000329490 0.21378639844
## [21] 0.14111368237 0.23186297483 0.29601755581 0.29336485560
## [25] 0.36276560985 0.53421634167 0.41002618666 0.30536282955
## [29] 0.75406967478 0.57026132394 0.70299558387 1.11431934970
## [33] -2.14164687950 0.22831008994 0.53968654931 0.52913420631
## [37] -0.53846008443 0.44463817520 0.30271805641 0.77162464170
## [41] -0.09091985909# standard error of estimator
# insoluable--bootstrapping may be the best solution here0.0.9 Spatial Point Pattern Analysis
#R code to test complete spatial randomness
#Source: Chapter 7 of Bivand, Pebesma, & Gomez-Rubio 2008
#Usage of this code requires citation of the original source
#clean up and set directory
rm(list=ls())
#load libraries
library(spatstat)
library(maptools)
library(rgdal)
library(maps)
#load three data sets
data(cells)
data(japanesepines)
data(redwood)
###Fun with Plotting###
#redefine as SpatialPoints objects
spcells<-as(cells,"SpatialPoints"); summary(spcells)## Object of class SpatialPoints
## Coordinates:
## min max
## [1,] 0 1
## [2,] 0 1
## Is projected: NA
## proj4string : [NA]
## Number of points: 42spjpines<-as(japanesepines,"SpatialPoints"); summary(spjpines)## Object of class SpatialPoints
## Coordinates:
## min max
## [1,] 0 5.7
## [2,] 0 5.7
## Is projected: NA
## proj4string : [NA]
## Number of points: 65spredwood<-as(redwood,"SpatialPoints"); summary(spredwood)## Object of class SpatialPoints
## Coordinates:
## min max
## [1,] 0 1
## [2,] -1 0
## Is projected: NA
## proj4string : [NA]
## Number of points: 62#convert to a unit square
spcells1<-elide(spcells,scale=TRUE,unitsq=TRUE); summary(spcells1)## Object of class SpatialPoints
## Coordinates:
## min max
## [1,] 0 1
## [2,] 0 1
## Is projected: NA
## proj4string : [NA]
## Number of points: 42spjpines1<-elide(spjpines,scale=TRUE,unitsq=TRUE); summary(spjpines1)## Object of class SpatialPoints
## Coordinates:
## min max
## [1,] 0 1
## [2,] 0 1
## Is projected: NA
## proj4string : [NA]
## Number of points: 65spredwood1<-elide(spredwood,shift=c(0,1),scale=TRUE,unitsq=TRUE); summary(spredwood1)## Object of class SpatialPoints
## Coordinates:
## min max
## [1,] 0 1
## [2,] 0 1
## Is projected: NA
## proj4string : [NA]
## Number of points: 62#plot the three data sets from Figure 7.1 of Bivand, Pebesma, and Gomez-Rubio 2008
pdf('cells.pdf');plot(spcells1); axis(1); axis(2); box();dev.off()## quartz_off_screen
## 2pdf('pines.pdf');plot(spjpines1); axis(1); axis(2); box();dev.off()## quartz_off_screen
## 2pdf('redwood.pdf');plot(spredwood1); axis(1); axis(2); box();dev.off()## quartz_off_screen
## 2###testing complete spatial randomness (Note:the code in the "plotting" section is unnecessary for this)###
#For the Japanese Pines
envjap<-spatstat::envelope(japanesepines,fun=Gest,nrank=2,nsim=99) #gets confused relative to "boot" library's version## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.#plot(envjap)
plot(y=envjap$obs,x=envjap$theo,type='l')
lines(y=envjap$lo,x=envjap$theo,lty=3,col='red')
lines(y=envjap$hi,x=envjap$theo,lty=3,col='red')
##Use the F function: distance from a point to the nearest event##
Fenvjap<-spatstat::envelope(japanesepines,fun=Fest,nrank=2,nsim=99)## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.plot(y=Fenvjap$obs,x=Fenvjap$theo,type='l')
lines(y=Fenvjap$lo,x=Fenvjap$theo,lty=3,col='red')
lines(y=Fenvjap$hi,x=Fenvjap$theo,lty=3,col='red')
###Redo for Redwoods###
envred<-spatstat::envelope(redwood,fun=Gest,nrank=2,nsim=99)## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.plot(y=envred$obs,x=envred$theo,type='l')
lines(y=envred$lo,x=envred$theo,lty=3,col='red')
lines(y=envred$hi,x=envred$theo,lty=3,col='red')
Fenvred<-spatstat::envelope(redwood,fun=Fest,nrank=2,nsim=99)## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.plot(y=Fenvred$obs,x=Fenvred$theo,type='l')
lines(y=Fenvred$lo,x=Fenvred$theo,lty=3,col='red')
lines(y=Fenvred$hi,x=Fenvred$theo,lty=3,col='red')
###Redo for Cells###
envcell<-spatstat::envelope(cells,fun=Gest,nrank=2,nsim=99)## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.plot(y=envcell$obs,x=envcell$theo,type='l')
lines(y=envcell$lo,x=envcell$theo,lty=3,col='red')
lines(y=envcell$hi,x=envcell$theo,lty=3,col='red')
Fenvcell<-spatstat::envelope(cells,fun=Fest,nrank=2,nsim=99)## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.plot(y=Fenvcell$obs,x=Fenvcell$theo,type='l')
lines(y=Fenvcell$lo,x=Fenvcell$theo,lty=3,col='red')
lines(y=Fenvcell$hi,x=Fenvcell$theo,lty=3,col='red')
###############################################################
#Source: Monogan, Konisky, & Woods (N.d.)
#Usage of the code or data requires citation of the original source
############## ESTIMATING MODELS #########################
#clean up and set directory
rm(list=ls())
options(scipen=10,digits=10)
#load libraries
library(spatstat)
library(maptools)
library(rgdal)
library(maps)
library(foreign)
library(mgcv)
library(rgeos)
#library(xtable)
library(akima) #for drawing image/surface plots
library(scatterplot3d)
#library(spdep)
checkCRSArgs <- function(uprojargs) {
.Call("checkCRSArgs", uprojargs, PACKAGE="rgdal")
}
#Load data
full<-read.csv('mnAir.csv')
#GAM model
gam.air<-gam(air~distances+s(longitude,latitude), data=full, family=binomial(link='logit')); summary(gam.air)##
## Family: binomial
## Link function: logit
##
## Formula:
## air ~ distances + s(longitude, latitude)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.0941008720 0.3020727202 3.62198 0.00029236 ***
## distances -0.0021650171 0.0009046893 -2.39311 0.01670643 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(longitude,latitude) 4.3121 5.969084 9.31141 0.15312
##
## R-sq.(adj) = 0.0318 Deviance explained = 3.11%
## UBRE = 0.32832 Scale est. = 1 n = 548for.comparison<-gam(air~distances+s(eastings,northings), data=full, family=binomial(link='logit')); summary(for.comparison)##
## Family: binomial
## Link function: logit
##
## Formula:
## air ~ distances + s(eastings, northings)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.100839924 0.301218429 3.65462 0.00025756 ***
## distances -0.002185183 0.000901764 -2.42323 0.01538311 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(eastings,northings) 4.435246 6.112045 9.52867 0.15314
##
## R-sq.(adj) = 0.0325 Deviance explained = 3.18%
## UBRE = 0.32787 Scale est. = 1 n = 548#draw figure of smoothed intercept
mn<-map("state","minnesota",resolution=0)
vis.gam(gam.air,view=c("longitude","latitude"),color='gray',
plot.type='contour',type='link',cond=list(distance=0),main="",
xlab="Longitude", ylab="Latitude",ylim=mn$range[3:4],xlim=mn$range[1:2])
lines(mn,col='red',lwd=2)
points(x=full$longitude[full$air==1], y=full$latitude[full$air==1], pch='+',col='white')
points(x=full$longitude[full$air==0], y=full$latitude[full$air==0], pch=1,col='white')
############# Testing for complete spatial randomness. If it doesn't fail, something is bizarre. ###################
#Draw a map of the 50 states
state<-data(stateMapEnv)
us.states<-map("state", fill=TRUE, plot=TRUE, resolution=0) 
IDs <- sapply(strsplit(us.states$names, ":"), function(x) x[1])
us.map.0<-map2SpatialPolygons(us.states, IDs=IDs,
proj4string=CRS("+proj=longlat +datum=NAD27"))
us.map<-spTransform(us.map.0, CRS("+proj=longlat +datum=NAD83"))
#extract minnesota
mn<-SpatialPolygons(Srl=list(us.map@polygons[[22]]),
proj4string=CRS("+proj=longlat +datum=NAD83"))
reach<-mn@bbox[1,2]
#prepare for CSR test
mn.trim<-mn@polygons[[1]]@Polygons[[1]]@coords
mn.coords<-dim(mn.trim)[1]-1
mn.bound<-owin(xrange=c(min(mn.trim[,1]),max(mn.trim[,1])),
yrange=c(min(mn.trim[,2]),max(mn.trim[,2])),
poly=list(x=rev(mn.trim[1:mn.coords,1]),
y=rev(mn.trim[1:mn.coords,2])))
mn.ppp<-ppp(x=full$longitude, y=full$latitude, window=mn.bound)## Warning: 9 points were rejected as lying outside the specified window## Warning: data contain duplicated points# There are 548 observations, but 8 are tossed for being illegal.
#look at it
plot(mn.ppp)## Warning in plot.ppp(mn.ppp): 9 illegal points also plotted
#conduct CSR tests
env.mn<-spatstat::envelope(mn.ppp,fun=Gest,nrank=2,nsim=99); env.mn## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.## Pointwise critical envelopes for G(r)
## and observed value for 'mn.ppp'
## Edge correction: "km"
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Significance level of pointwise Monte Carlo test: 4/100 = 0.04
## .....................................................................
## Math.label Description
## r r distance argument r
## obs hat(G)[obs](r) observed value of G(r) for data pattern
## theo G[theo](r) theoretical value of G(r) for CSR
## lo hat(G)[lo](r) lower pointwise envelope of G(r) from simulations
## hi hat(G)[hi](r) upper pointwise envelope of G(r) from simulations
## .....................................................................
## Default plot formula: .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 0.15055]
## Available range of argument r: [0, 0.41665]#env.mn<-spatstat::envelope(mn.ppp,fun=Gest,nrank=2,nsim=99,savepatterns=TRUE); env.mn
plot(y=env.mn$obs,x=env.mn$theo,type='l',xlim=c(0,1),ylim=c(0,1))
lines(y=env.mn$lo,x=env.mn$theo,lty=3,col='red')
lines(y=env.mn$hi,x=env.mn$theo,lty=3,col='red')
Fenv.mn<-spatstat::envelope(mn.ppp,fun=Fest,nrank=2,nsim=99); Fenv.mn## Generating 99 simulations of CSR ...
## 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38,
## 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76,
## 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99.
##
## Done.## Pointwise critical envelopes for F(r)
## and observed value for 'mn.ppp'
## Edge correction: "km"
## Obtained from 99 simulations of CSR
## Alternative: two.sided
## Significance level of pointwise Monte Carlo test: 4/100 = 0.04
## .....................................................................
## Math.label Description
## r r distance argument r
## obs hat(F)[obs](r) observed value of F(r) for data pattern
## theo F[theo](r) theoretical value of F(r) for CSR
## lo hat(F)[lo](r) lower pointwise envelope of F(r) from simulations
## hi hat(F)[hi](r) upper pointwise envelope of F(r) from simulations
## .....................................................................
## Default plot formula: .~r
## where "." stands for 'obs', 'theo', 'hi', 'lo'
## Columns 'lo' and 'hi' will be plotted as shading (by default)
## Recommended range of argument r: [0, 0.42564]
## Available range of argument r: [0, 0.42564]plot(y=Fenv.mn$obs,x=Fenv.mn$theo,type='l',xlim=c(0,1),ylim=c(0,1))
lines(y=Fenv.mn$lo,x=Fenv.mn$theo,lty=3,col='red')
lines(y=Fenv.mn$hi,x=Fenv.mn$theo,lty=3,col='red')
#Very clear patterns of clustering on both counts.
#Fit a point process model
fit.mn<-ppm(mn.ppp,~1); summary(fit.mn)## Point process model
## Fitting method: maximum likelihood
## Model was fitted analytically
## Call:
## ppm.ppp(Q = mn.ppp, trend = ~1)
## Edge correction: "border"
## [border correction distance r = 0 ]
## ---------------------------------------------------------------------------
## Quadrature scheme (Berman-Turner) = data + dummy + weights
##
## Data pattern:
## Planar point pattern: 539 points
## Average intensity 21.1 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 353 vertices
## enclosing rectangle: [-97.2252121, -89.47309113] x [43.49322891,
## 49.38323212] units
## Window area = 25.5331 square units
## Fraction of frame area: 0.559
##
## Dummy quadrature points:
## 64 x 64 grid of dummy points, plus 4 corner points
## dummy spacing: 0.12112689018 x 0.09203130007 units
##
## Original dummy parameters: =
## Planar point pattern: 2367 points
## Average intensity 92.7 points per square unit
## Window: polygonal boundary
## single connected closed polygon with 353 vertices
## enclosing rectangle: [-97.2252121, -89.47309113] x [43.49322891,
## 49.38323212] units
## Window area = 25.5331 square units
## Fraction of frame area: 0.559
## Quadrature weights:
## (counting weights based on 64 x 64 array of rectangular tiles)
## All weights:
## range: [0.000429, 0.0111] total: 25.5
## Weights on data points:
## range: [0.000429, 0.00557] total: 1.47
## Weights on dummy points:
## range: [0.000429, 0.0111] total: 24
## ---------------------------------------------------------------------------
## FITTED MODEL:
##
## Stationary Poisson process
##
## ---- Intensity: ----
##
##
## Uniform intensity:
## [1] 21.10985456
##
## Estimate S.E. CI95.lo CI95.hi Ztest
## log(lambda) 3.049739972 0.04307304923 2.965318347 3.134161597 ***
## Zval
## log(lambda) 70.80390237
##
## ----------- gory details -----
##
## Fitted regular parameters (theta):
## log(lambda)
## 3.049739972
##
## Fitted exp(theta):
## log(lambda)
## 21.10985456plot(simulate(fit.mn,2))## Generating 2 simulated patterns ...1, 2.
#beautiful, the simulation (also used in "envelope") respects the polygon boundary